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ABSTRACT 


This work is concerned with two inter-related problems: (1) the 
theory and application of algorithms for the discrete integration of the 
acoustic wave equation in large, multi-dimensional, variable-parameter 
domains and, (2) the automatic synthesis of adaptive acoustic arrays 
for the detection and estimation of received acoustic signals which are 
incident from such domains, 

It is shown that the steady-state integration of large domains may 
be partitioned into subdomain integrations for slowly-varying variable- 
parameter domains. A new method of integration based upon the Fast 
Fourier Transform is given, A new method for obtaining closed form 
coefficient expressions for the Fast Fourier Transform is shown and 
illustrated, 

The results of a general computer-display simulation program for the 
Widrow feed-forward algorithm are given. Several limitations and 


possible modifications to the Widrow procedure are given. The use of 


Kalman-type filters as an alternative method is introduced. 
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I. INTRODUCTION 

The detection, estimation, and distribution of acoustic energy signals 
propagating in inhomogeneous seawater are subjects of continuing interest. 
These subjects reach a focus of inter-relation in the design of efficient 
acoustic receiving arrays; such arrays Should match the propagation 
solution properties of the location where the receiver is to be placed. 

Past investigations into the trajectory and strength of acoustic 
Signals in seawater have relied almost totally upon minimum-time plane- 
wave ray-path techniques or approximate normal mode solutions to crude 
stratification models. The theoretical and computational economy of 
these methods is their main asset. Rarely is a single model sufficient 
to represent faithfully the propagation characteristics of long acoustic 
paths involving multi-dimensional energy spreading or inhomogeneous 
medium phenomena such as convergence Zones. 

A more fundamental approach to such problems is the direct inte- 
gration of the appropriate variable-coefficient wave partial-differential 
equation with any receiving array present being treated as part of boundary 
constraints. Unfortunately, no general theoretical methods exist for 
treating variable-coefficient partial-differential equations. However, for- 
seeable order of magnitude increases in digital computation speed and 
memory capacity make the direct finite-difference numerical integration of 


the wave equation a distinct possibility. 
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This paper attempts to make contributions to both the numerical 
integration techniques relative to solution of the large-domain variable- 
parameter acoustic wave equation and the theory and practice of adaptive 
acoustic antenna arrays. 

In chapter II the acoustic wave partial-differential equation with 
variable coefficients is derived in properly posed form from first prin- 
ciples assuming infinitesimal deformations and that the coefficients are 
Slowly varying and time independent. The appropriate boundary condi- 
tions are also discussed with particular emphasis being devoted to the 
matched boundary condition. 

In chapter III methods for finite difference integration of discrete 
models for partial-differential equations are introduced. The problems of 
instability and convergence are illustrated for the one-dimensional 
diffusion equation which is the simplest partial-differential equation 
of common interest which has been extensively studied, 

Continuous-discrete, explicit, and implicit solution algorithms are 
then studied and the results of actual computer calculations relative to 
required computation time, accuracy and stability are presented and dis- 
cussed, Along with these results, a new method for semi-implicit 
integration of the three-dimensional wave equation using the Fast Fourier 
Transform algorithm of Cooley and Tukey is given. The method appears 
faster than all other methods tested with the exception of the successive 


approximation tri-diagonal method of Douglas and Gunn. 
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The magnitudes of the spacial increments Ax, Ay, Az for proper 
wave equation simulations are obviously highly important; it is shown that 
the highest desired propagation frequency places a maximum value on 
these increments which is considerably less than the appropriate Shannon 
restriction. 

The problem of integrating much larger spacial domains than could 
possibly be contained in any reasonable computer memory is taken up in 
chapter IV. General transient problem integrations appear to be very 
difficult to treat in such cases; the breaking up of the large integration 
domain into smaller memory requirement domains which can be integrated 
is frustrated in the variable parameter case by the question of which small 
domain is to be integrated first. It is shown however, that the steady 
state solution problem can be treated in such a manner for the slowly- 
varying parameter case by making use of the matched impedance boundary 
condition of ordinary transmission line theory. Since the problem of calcu- 
lating steady state acoustic propagation loss is of this type the result is 
Significant. The chapter concludes by presenting the results of several 
computer wave-equation integration problems. 

The Fast Fourier Transform algorithm is a very efficient digital computer 
algorithm for computing the discrete Fourier transform of a finite-length 
sampled input function. Although a few discrete finite transforms can be 
obtained in closed form by use of the geometric sum formula, no general 
techniques for obtaining the explicit closed-form coefficients presently are 
known. In chapter V a new technique for obtaining such forms is given 


and two examples are worked as illustrations. The method is based on 


7 


taking first-order backwards finite differences and is quite general, 

The subject of acoustic detection and estimation arrays is taken up 
in chapter VI. From the aspect that the local medium propagation proper- 
ties may vary with the time and as a consequence the most efficient 
detector should also, the concept of automatic self-adapting arrays appear 
especially appealing. Therefore, after initially considering arrays from the 
classical optimum detection point of view, adaptive arrays are studied. 
A quite sophisticated general purpose interactive-display digital 
computer program has been used to study several possible automatic 
synthesis adaptive antenna arrays which may be placed on the boundary 
of a large wave-equation integration domain. General results are 
discussed and particular emphasis is devoted to possible modifications of 


existing adaption algorithms in order to increase their adaption speed, 
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II, THE PROPAGATION EQUATIONS 

A. THE UNBOUNDED MEDIUM EQUATIONS 

The problem of propagation of acoustic waves in salt water will be 
assumed to satisfy an acoustic wave equation; however, the solution for 
domains with variable parameters in more than one dimension is very 
complicated in general so that it is advisable to begin the analysis with 
a derivation of the medium equations from first principals in order to make 
clear the assumptions involved. 

For infinitesimal acoustic disturbances in fluids the volume bulk 


[1] 


modulus x is related to the sound speed c by the relation 
x = pc (1) 
where p is the density of the fluid. For the range of conditions dealt 
pau nere, the density 9 will be considered #6 be a censtant. The acoustic 
(2 |] 


speed c, however, is a function of both temperature and pressure. 
Notwithstanding these complications, it is usual to assume that the 
resulting equation of state 

C Tem pyr) (2) 


can be reduced to a parametric, Single-valued function of position 


Care sos, 72) (3) 
= ¢ (x) 
which is constant with respect to time for periods long compared with the 
eee (2] LJ | 
periods of oscillations of interest, By Hookes law, the compression 


of particles relative to an increment change in pressure, p(x,y,zZ,t) is 


p (x,y,2,t) =-x 7° a (x,y,2,t) (4) 


ig 


where d(x,t) is the vector relative displacement of the particles. Taking 
the partial derivative of (4) with respect to time gives 


o 
at 


= -x7-v(&,t) (5) 
where v(x,t) is the vector acoustic disturbance velocity. The two depen- 
dent variables p,v are further related by Newton's second law which relates 


the vector acceleration of a sample mass volume to the pressure gradient 


d 
op = - 7p. (6) 


It is usual to further simplify (6) by assuming that in the total veloe= 


ity derivative 


a en a Amn ae 2 (7) 
the second term is of second order for all but the very highest possible 
frequencies so that 

Ov 
mi ao (8) 
The two equations (5) and (8) constitute the fundamental first-order 
partial-differential equations to be solved in the absence of dissipative 
losses. The effects of dissipative losses usually are included empiri- 
cally after solution of the lossless equations (5) and (8). ui] 
Equation (8) is separable into three scalar partial-differential 


equations so that with (5) there are a total of four first-order partial 


differential equations to be solved: 
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ap Ov Ov Ov 


—— ey —_z (9) 
ot ies Ox 3 Oy i 82 
Ov 1 ep 

5 
a ] 
ot Pp ox (10) 
Ov lor 
—%=-— — (11) 
ot p 9y 
Mg en = is : (12) 
ot eee 


The bulk modulus x is a function of position since by (3) c is a function 
of position and p has been assumed constant. The set of equations (9,10, 
11,12) are a mixed set of equations in that they are functions of two 
dependent variables p and v. They may be unmixed by taking the partial 


of (9) with respect tot, 








a” p i —- 
sis So oe O ek 
8 
en eg ee p 
0 
ee, vos (14) 


This is a scalar wave equation for the pressure, The velocity v 
satisfies an associated type of equation; if the partial of (8) with respect 


to tis taken and continuity is again assumed, 





o“v 1 
st2 ad (= i 
i 
=- = 7(-x7°Vv) 
fe 
= U(Tev) 7m. (15) 


Z1 


The problems of interest here will be confined to unbounded fluids which 
are assumed unable to support shear stress waves; such being the case 
the curl of v must vanish and hence by the vector identity 

7(7-v) = 7x (vxv) + Vv (16) 
(15) becomes 

o° Vv x l 

Se ad (JVev) TH, (17) 
The second term in (17) is usually neglected as being of second order 
since 7* v is assumed small for the reason previously stated and in 
addition 7x can be made also small by proper space partitioning in 
computer analyses, In ray-path solutions the second term can become 
important in certain cases, i] 

As previously stated, the equations (9,10,11,12) are fundamental; 

the derived set (14) and (17) are no less valid but their domain of solutions 


[4] 


is larger. It is possible for certain pathological solutions to exist 
satisfying (14) and (17) under a given set of boundary conditions but 

these solutions none-the-less will not satisfy (9,10,11,12). Thus only 
(9,10,11,12) are both necessary and sufficient, The pathological solutions 
to (14) and (17) will almost never occur in practice if (14) and (17) are 
solved analytically since the investigator will invariably pick character- 
istic forms in separation procedures which are physically regular in behavior; 
the problem is not academic however in computer solutions of the equations 
Since poorly defined (real world) boundary conditions could possibly lead 


to successive approximate solutions converging to these pathological 


solutions. For all of these reasons the integration techniques in this work 
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will apply to the equations (9,10,11,12). 
The set of equations (9,10,11,12) is also properly posed, i.e., a norm 
can be defined in an appropriate Banach space which is bounded, 4] 
In implicit algorithm solutions to (9,10,11,12), itis very con- 
venient to introduce new variables to make the equations (5) and (8) 
completely symmetric. Let 
aa (18) 


p= Bp’ (19) 


then to first order 


a = C7 (20) 
av" 
at == CVJe p (21) 
where 
ae (22) 
= constant 


l 
a (23) 
a(x) oc (x) 
The first order assumption in (20) and (21) assumes that a (x) is a slowly 
varying function of position. In addition to symmetry the new variables 


have a much reduced dynamic range which considerably increases the 


accuracy of computer computations. 


B. BOUNDARY CONDITIONS 

The equations (5) and (8) require a complete, consistent set of bound- 
ary conditions in order to guarantee that their solution will be unique and 
properly posed. The usual discrete methods of solution, from a practical 


if not a mathematical point of view, almost always are limited to bounded 
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domains. This excludes several important practical problems of interest 
including wave propagation in semi-infinite media, For the special case 
in which the steady-state propagation properties are desired it will be 
shown that by an appropriate spacial partitioning of the volume domain 
of interest together with an extension of the matched boundary condition 
of one-dimensional transmission line theory this problem may be solved, 

The boundary condition at a fixed wall requires that the normal com- 
ponent of velocity vanish 

n * v (bnd,t) = 0 (24) 

where nis a unit normal to the boundary surface. 
An equivalent, not independent, condition on the pressure is from (6) 


poe 2 fh donee (25) 


At a free surface, if surface tension forces and atmospheric pressure may 
be ignored, the pressure must vanish 
p (surface, t)=0. (26) 
The matched boundary condition relates the pressure and velocity 
at a point and is an example of a mixed boundary condition. The condi- 
tion is most fundamentally applied to one dimensional problems and is 
precisely reflectionless for constant parameter applications. The condi- 
tion is 
rey 7 0% 27 
where |v| represents the magnitude of the velocity. The condition can be 
approximately extended to variable-parameter media as well as multiple 


dimensions. 
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C, COORDINATE SYSTEMS 

The equations (5) and (8) are in vectorial form even though derived 
from rectangular coordinate considerations; they hence hold true for any 
orthogonal curvilinear coordinate system. The choice of coordinate 
system depends upon the type of problems to be investigated as well as 
the domain size and approximations to be utilized. For example, physical 
problems inherently possess propagation properties which are three- 
dimensional, but frequency, source location, parameter variations and 
fixed boundaries may allow the propagation process to be confined to 
only one or two dimensions. 

The discussions in this work are confined to rectangular and 
cylindrical coordinates but the technique often can be extended ina 
simple manner to spherical coordinates. The theory is usually illustrated 
for simplicity in one dimension but the computer simulations that have 


been carried out are always at least two-dimensional, 


D. THEORETICAL SOLUTIONS TO THE PROPAGATION EQUATIONS 

To measure the accuracy of computer solutions to (5-8), theoretical 
solutions are desirable. In those cases where the parameters are 
constants, the ordinary theoretical acoustic solutions to the wave equation 
in homogeneous media suffice; in the general case of inhomogeneous media 
exact theoretical solutions are almost non-existent except for some very 
special cases. The case of linear periodic-parameter media has received 

[5 ] 


extensive treatment and constitutes the most useful theoretical solution 


known as a comparison standard for finite discrete models of the wave 
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equation. In chapter JII it is shown that such models are inherently 
dispersive, whereas the true wave equation in unbounded media is perfectly 
non-dispersive. However, if the partitioning is sufficiently small, the 

low frequency or Rayleigh limit approaches precisely the non-dispersive 
Linnie 

The step function is used as a test function in this work to verify 
that the Rayleigh limiting velocity is approached and that dispersive wave- 
form distortion is not prohibitive. 

Finite difference approximations to the wave equation are band- 
limited inherently so that the infinite slope at the discontinuity of a step 
function will never be represented exactly in computer solutions. The 
maximum attainable slope of the solution is most easily found by use of 


[6 | 


the Fourier transform, Assume that the resulting solution at some 
arbitrary point a and time t is band-limited to | f|< B hertz. Then by 
the Fourier transform the resulting pressure is represented by 


B 
p@t=f, G@f 


where G(a, f) is the spectral resolution function. Hence 


ee (28) 


a <2mB + Max (|G(a,f)|) 2B. (29) 
Assuming that the pressure wave is an outward-going wave in the x 
direction only then the velocity in the x direction, by (27) is directly 
proportional to the pressure and hence by (9) 


ap QV 


ok 
at ti(‘<é‘ SC (30) 
= cee Os 
JV%0 Ox 


or 
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ej = 2 |B 


< fe. 4a B° Max (lG(a,£)| Nee 


Thus the frequency band-limitedness of discrete models of the wave 


(31) 


equation has as a direct consequence the result that the spacial repre- 
sentation of waveforms is also spacial-slope limited. In essence, then, 
transient-solution curves for step-function time input to discrete models 

of the wave equation will exhibit finite-maximum spacial-slope properties 
as the solution is observed at various fixed times, This has been observed 
in computer solutions and will be discussed in the chapter on discrete 


models for the wave equation, 


Zu, 


III, DISCRETE INTEGRATION 

A, THE DIFFUSION EQUATION 

The purpose of this chapter is to introduce discrete integration models 
for partial differential equations Suitable for programming on digital, 
analog, or hybrid computing machines. The relevant theory of such models 
is extensive. The heat diffusion equation in one space dimension is the 
simplest partial-differential equation possessing known theoretical 
solutions and which also illustrates the problems of instability and 
convergence common to discrete difference methods of solution. The wave 
equation, being a hyperbolic partial differential equation of the second 
order, is inherently more difficult to solve than the diffusion equation 
which is parabolic and of first order so that the relevant theory will be 
illustrated for the diffusion equation with necessary modifications for 
the wave equation being taken up in section C. 

I theoretical solution 

The diffusion of heat through a one-dimensional solid of thermal 
conductivity n and heat capacity a is given by the parabolic partial 


differential equation 


ag _2 7, at 
At ax A ee ie) 


where T(x,t) is the temperature of the solid at position x and timet. 
Both the heat capacity and conductivity can in general be functions of 
position and time; it shall be assumed here that they are constants so 


that (32) becomes 
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aT a°T 


—_ = —> 3 
at 2 ax” (33) 
where 


aie 
a 


a (34) 


For a domain assume the real line osx <7 and thatt 20 


Further assume the given boundary values 

T (x,0o) = g(x) (35) 
and 

T (o,t) =T (7, t)=o0. (36) 


As an example suppose 


p(x) = cx , DEE 
= c(7- x), a <x <7 , (37) 


The theoretical solution to this boundary value problem is obtainable 


in convergent infinite Fourier series form; the assumed solution 
_ 2 jmx 
T (x,t) = yA. (tHe (38) 
upon substitution in (33) requires that 
dA 


ye a - o (im? aA ber =o (39) 


which can only be true for arbitrary values of x provided 


dA. 
eo 2 = 
a + m oA a (40) 


The solution to (40) which is bounded for allt2o is 
2 
- -om*t 
Aft) = A_o)e ; (41) 


Hence (38) may be written as 


jmx-m* ot 


T(x,t) = B__ Am (o)e (42) 


Zo 


The coefficients A, (o) are obtained by the usual coefficient formula from 


the distribution of temperature at t = 0: 


Sei -jmx 
A_(o) =|" ~ (x)e dx 
=o, m even 
= ope ; m odd, (43) 


This theoretical solution has been evaluated for different values of 
the time by summing the series (42) with coefficients (43) on an IBM 
360/67 to six decimal places; the solution is shown in Fig. l. 

The use of the Fourier series method of solution is limited strictly to 
linear constant-coefficient partial-differential equations whereas the 
finite difference methods to follow can be applied to much more general 
types of problems; the theoretical solution however does give a standard 
of comparison for any proposed finite-difference integration scheme 
when such a scheme is applied to linear constant-coefficient problems, 

2. Discrete-Difference Equation - Instability 

A great number of possible discrete models can be proposed to 
integrate numerically the equation (33) with the domain of (35,36). The 
majority of methods discretize x and t into some kind of grid space in 
order to convert the solution of the partial-differential equations or a set 
of algebraic equations. The grid may be uniform or non-uniform in the 
coordinates mA&x, nAt; the crucial questions to be answered for any 
scheme, however, usually are as follows: 


A. Is the integration scheme stable as the number of points in 
the mAx, nAt solution space become infinite? 


B. Does successive refinement of the solution by increasing the 
number of grid points in the solution space produce a solution 
which converges uniformly to the theoretical solution? 
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C. From a practical point of view the number of grid points must 
always remain finite, At what rate as Ax, At — o does the 
error between the finite difference approximation and the exact 
solution tend to zero? 

The answers to these questions are exceedingly difficult in 
general; usually each case must be treated separately, especially in 
variable coefficient or non-linear problems. The known theory relative to 
questions A and B is fairly complete, especially for linear problems; the 
present status of knowledge about question C is quite inadequate since 
relative few exact solutions are known for variable coefficient problems, 

The simplest and most commonly applied discrete-approximation 
scheme to the problem solved above is to discretize the space variable 


x into M sub-intervals Ax where 


Lx =— , (44) 


m fight One n nt n 
At ~ At? es mT bas 


mos Ie. we, M-l, n=0,1,2,.22 


The boundary conditions given in (36) become 


i — Tlo,n At) 
O 
=o n=o,l (46) 
n 
T.. = T(MAx, nAt) 
M 
= Oo n=o,l (47) 
and the initial cordition (37) prescribes 
T° = T(mAx,o) 
m 
= ¢~ (m Ax) Onsen ie (48) 
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Nothing has been said about the magnitude of the time increment At; 
this is the crux of the problem as is indicated in the plots of the computer 
solution to (45) shown in Figs. 1, 2(a,b,c). For small values of At the 
solution predicted by (43) agrees quite well with the theoretical solution 
shown in Fig. 1; when the value of At becomes large the computed 
solution becomes unstable as shown in Fig. 2(a,b,c). This is not caused 
by roundoff or truncation errors of the digital computer (double precision 
calculations have no effect) but is a direct result of the instabilities 
inherent in (45). 

3. Exact Solution of the Difference Equations 


[36] 


Von Neumann was the first investigator to successfully explain 
the instability inherent in the system (45) by means of Fourier analysis. 


To apply the method to (45) assume a series of the form 


7: - > Ber eS ee (49) 
mM k= -o@ 
is substituted in (45). After some algebra it can be seen that (49) 
satisfies the system provided 
Zo At 
T(k) =1- (Reds (1 - cos (kK Ax)) . (50) 


Further, if the coefficients B, are set equal to the corresponding A. (o) in 
(43), the solution (49) is exactly the same as the previous Fourier solu- 
tion (42) except for the time dependence function. In (49) only the values 
at the actual mesh points mAx, nAt of course have meaning. Assuming 
then, that both the Fourier solution (42) and the solution 


is = > A, (0) T” _)k(m 4x) (51) 
k= —- © 
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are to be evaluated at the mesh points mAx, m=o0,1,..,M, the essential 
difference in the solutions occurs in the amplitude of the time factors 


“0 m*t . The factors represent the time dependent perfor- 


n 
(rT (k)) and e 
mance of the individual Fourier harmonics in respective solutions and 


because of linearity are independent of each other. The exponential factor 


z: 2 
3 opt clearly decays for increasing time whereas the factor 
n 2 l 4 2 Ph 
(rt (k) ) = (I-k° oAt + 95 kK oAt(Ax)" - ... ) (52) 
tends to zero for increasing n only if 
Pig) = 1 (53) 


If for any harmonic the inequality (53) is not satisfied, that particular 
harmonic will grow as time increases leading to the instability shown in 
ice 2 (a, D,C) . 

From (50) t is real for all values of k and never more positive 
that + 1. The only possibility for instability thus must result from those 
cases where fT is less than -1. The worst condition occurs for the 
particular value of k = kK. such that 

cos (k Ax) = -l, (54) 
which for stability restricts 


ZOuNt 
(Ax)* 


Conceptually, the equality in (54) occurs for those harmonics of 


<a liens (S'5) 


wave-number Ko whose wave-lengths are of the order of magnitude of 
twice the discrete spacial grid spacing Ax. If these harmonics have 
amplitude exceeding the inequality restriction in (52) they will become 
unacceptably large as the time progresses, It has been shown that the 
inequality (55) is both necessary and sufficient for stability of this 


discrete integration scheme, 33 


4. Refinement of Mesh - Convergence 


If the domain of x does not change, itis clear that the magnitude 
of Ax will decrease as M increases, i.e., as the grid or meshis 
refined. Question B previously stated in essence asks whether, as the 


grid is refined the limit 


n 


a | T (mAx, nAt) - Ta, 


ve ao (56) 


holds true. In (56) T(mAx, nAt) is the value of the temperature 
calculated by substituting the values mAx, nAt in the exact theoretical 


[37 ] 


solution of the diffusion equation. Hildebrand has proven that the 
restriction (55) is sufficient to guarantee the limit (56). Further it can be 
stated that if a like analysis can be performed on other discretization 
algorithms for solution of the diffusion equation which lead to analogous 
expressions equivalent to (55) the sufficiency of convergence follows, 
So eRales iofeConvergence 

The problem of rate of convergence (Question C), usually must be 
handled on an individual basis and experience in actual computation is 
still the fundamental requirement. 

In the case of the wave equation it will be shown in Chapter III 
that the minimum number of spacial increments required is set by the 


upper frequency limit for which the propagation properties are desired 


and by the characteristic acoustic velocity of the medium, 


B. IMPLICIT ALGORITHMS AND THEIR SOLUTION 
The restriction (55) on the relationship between At and Ax for 
stability is very severe; if Axis refined by a factor of *#, At must be 
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reduced to = of its previous value to guarantee stability. The solution 
algorithm (45) is called an explicit’ algorithm in that it explicitly 
predicts ‘aia directly. Implicit solution algorithms are also available 
which require the solution of a set of simultaneous equations in order to 
find oa . The more involved computations necessary are more than 
offset by the fact that the stability criterion may be much less severe, or 
in some cases, the system of equations may be unconditionally stable, 

In the latter cases the magnitude of the increment At is independent of 
Ax and can be set purely on the basis of the computing machine trunca- 
tion error involved or the fineness of the final solution grid desired, 

In general, if the partial differential equation is linear and only 
nearest-neighbor interactions are involved in the corresponding explicit 
algorithm, the choice of an implicit algorithm is indicated since the 
solution of the resulting linear equations involves inverting a nearly 
diagonal matrix, 

Reference 4 lists a great number of implicit schemes for the diffusion 


equation, some of which are unconditionally stable. Atypical example 


of an implicit scheme is the equation 


n+l n 
- 7 a Oo n+l n+l n+l 
At Wie iB veel : ae " TaD (57) 


n n n 
tte = + 
(1-8) oe ee Face | 
where f is an adjustable parameter. If B = o the scheme reduces to 
the explicit scheme considered previously; if f= 4a centered time-difference 


equation results which is implicit. The equation (57) can be investigated 


as previously for stability; the condition equivalent to (55) is 
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20 At il 





(xe * 1-28 ' if os B<% 
(58) 
no restriction it ee Pp ees 
[7] 


Hence, for example if 6=#% the equation is unconditionally stable. 


It can be written in the form 


n+l pinata laa meet 
- ny re At? Cel 
(5:3) 
gn ee n n 
= ones 2e)T + CT 41 + CT -l 
where 
— v0 Siete a 
C — 9 (Ax)* (60) 


The equation (59) contains M-1 unknowns, i.e., all of the variables 
on the left side of (59) minus the two known boundary conditions supplied 
by (46) and (47). The right side of (59) contains all known values so that 
the M-1 equations are sufficient for solution of the system by any 
of the standard methods applicable to the simultaneous solution of linear 
equations, 

The equations (59) can be written in the general form 


+ 
per’ = pene (61) 


in which A and B are tri-diagonal matrices. For the one-dimensional case 
there exists a particularly efficient adaption of the Gauss elimination 


algorithm for the solution of this system which will be discussed in section 


C on the wave equation, 


C. DISCRETE MODELS FOR THE WAVE EQUATION 
The wave equation, in component symmetrical form given by (20) and 
(21), is inherently more difficult to solve than the diffusion equation 


considered in section A. To begin with itis of second order in the time 
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as well as in space and hence requires twice as many boundary values and 
at least twice as much computation to solve, Further itis a hyperbolic 
equation and the solutions need not be analytic for all values of the 


9] 


independent variables as is the case always for the diffusion equation. 


[4] 


The formulation (20) and (21) constitute a properly posed system in that 


p’ and v’ constitute components of a general vector u which is dense and 


[8] 


possesses finite norm in the appropriate Banach space. 
1. Continuous - Discrete Integration Approach 


The simplest conceptual model for integrating the wave equation 
numerically is obtained by approximating the partial spacial derivatives 
occuring in (20) and (21) by finite differences. Although the approach is 
limited in digital applications by very severe restrictions on the time 
increment Atin order that stability be maintained, the approach is 
fundamentally valid and easy to apply being successfully used by many 
investigators. To be efficient computationally the procedure requires a 
digital computer having a great deal of parallel processing capability or 
an analog type computer. A general hybrid machine of large capacity 
would probably be ideal but as yet the amount of operational amplifiers 
and electronic switches required have prohibited the economic construction 
of such a machine. With the advent in the near future of integrated circuit 
technology this approach may become dominant, 

The result of replacing the partial spacial derivatives by finite 
difference approximations at a set of spacial grid points is to produce a 
set of ordinary-differential equations for the pressure and velocity at 


at these grid points, To illustrate the method assume in (20) and (21) 
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that the propagation is in the x direction with no variations inthe y and 


z directions. The one-dimensional first-order equations in normalized form 


are then 
Oo ciel Ni 
at = Ox (62) 
av 
ae 9p 
at “x (63) 


where the superscripts have been dropped. The equations are completely 

, ee ; Sees , 
analogous to the electrical transmission line equations in which 
pressure is analogous to signal voltage and velocity is analogous to 
signal current. (In the un-normalized form of (62) the coefficient c would 
be replaced by the reciprocal of the series inductance/length.) 


Conceptually the single variables p, im are replaced by 


column vectors v and p 


Po (t) 

p(t) = 
: (64) 
Py, 4) 
v(t) 
«O 

Eye (65) 
V4, (t) 

where 
P (t) = p(m Ax, t) (66) 
v_, {t) = v(m POL). (67) 


If the domain of x is of length £ and itis assumed subdivided into M 


subintervals then p and v are M + 1 dimensional and the system of 
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equations (62), (63) are replaced by the 2(M - 1) equations, 





1 Ymtl ~ Yn-1 

“ae =—- =- Cc (m Ax) 2 AX ’ m=1, oe « PRM 168) 
ae a= Cine Pm+l ee nl , M=l, ee, M-1 (69) 
dt 2X 


and four equations for Po Pa! Vas» and ve . These last variables are 
usually separated from (68) and (69) since their corresponding equation 
forms are dependent upon the appropriate boundary conditions, 

The boundary conditions appropriate to the system (68), (69) are 
obviously 2(M + 1) in number: the 2(M + 1) initial values if the system is 
completely homogeneous or 2(M - 1) initial conditions for the equation (68), 
(69) and pot), p(t), v(t), vas (t) if the system is assumed excited by 
its boundary values, The solution to this system is ideally carried out on 
an analog computer using separate interconnected integrators for each of 
the variables P Ce), va (t); the number of integrators is large however, 
and the setup and scaling problems involved are very cumbersome to 
handle for all but very small values of M. A hybrid machine with an 
automatic inter-connection matrix, potsetting and error control would be 
theoretically the best solution but no large machine meeting these require- 
ments presently exists or can be dedicated to so restricted a set of 
problems, Thus digital computers are usually employed to sequentially 
integrate such systems. The process is slow (relative to true problem time) 
even for the fastest, large-core machines now in existence, However, it 
is accurate and refinement of grid procedure can be made automatic. 


[9] 


Digitally, any of the common methods may be employed to 


integrate the ordinary differential equations (68), (69) but the accuracy is 
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clearly limited to the order of ( At+ Ax). The use of a small time incre- 
ment At is also required since stability requires that At/(Ax)* be bounded 
as Ax, At tend to zero. The simplicity of the scheme together with the 
fact that only nearest-neighbor interactions are involved make the scheme 
quite useful if real-time computations are of secondary interest. 


2. Explicit Scheme 


An explicit scheme which is slightly more accurate as well as 


stable (provided (c At/ Ax) <1) is 10] 

n+1 _ ot 

a mm ee es n ha n 
At «Ax va “m-3 ) (70) 

n+1 n 

m-2  m-3 
; See eee Pa 
At — A Do Die ; 


Here fractional subscripts (purely a notational device) have been intro- 
duced in order to avoid double space intervals 2 Ax. The scheme (70, 


71) can be arranged into purely explicit form as 


n+] n n n 
Pa mel tlt inte = Ringe! (72) 
n+ 1 n n n 
el = ae + C (D - P_-1) (73) 
3 n _ n Nn 
+C m=3 avy + Vin-3/2 ) 
where 
c At 
C = Ax ° (74) 


This scheme has been extensively used and can be generalized to three 


dimensions. The result is the system 
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p(jam,n,r+1) = p(t,.myngz) 
+t vy (2+ #,m,n,t)-v,(£-%,m,n,r) ) 
+ - (v (648 ;n,Y) - ya £,m-%,n,r) ) 


a S (v_( L /m,nt+s a) a a £,m,n-¥ /t) ) (75) 
-_ i = om 
veel L *,m 7H ie) we. | 4 + ,m,n,r) 
a Cy. (p( h Pisin, b) - p( £-1,m,n,r) ) 
2 st = = 
+v - (v( £-% ,m,n,r) av £- 1,m,n,r) 


+ v__ ues 3// eal ste) 


(76) 
a 3 = - 
WN £,mrg ,n,rtl1) Maal £,m-# ,n,r) 
a cee Mitr, ty, 5) p ( 4 ,»m-1,n,r) ) 
= 1 
+ v (£,m-% ,n,r) - 2v (£,m-1,n,r 
C2 (wv (4 pm-$ yn,t) ~ 2v | ) 
+ £,m-3/2 
Vek 7m 37 pa, t) ) 7) 


v (2% ,m,n-§ ,r+1) =v ( £ ,m,n-g ,1) 
+ (pC £,m,n,r) - p(£,m,n-1,r) ) 
2 
ta Gli (v_ £,m,n-% ,f) 2v_( £,m,n,-1,r) 
+ v_( £,m,n-3/2,r) ) 


(78) 


where 


p( 4opmynys) = splewAx amAyen Az ,nAtle 


ee — At _ cAt 


x Ax : My Ay : C z  4x£Az : (79) 


Computational experience on a CDC 6500 computer with this scheme 


has shown that: (1) The error is only slightly less than the scheme (67), 


(68); (2) The scheme is considerably more stable than (67), (68) 
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especially when the coefficients Ax, Ay, Az are not equal and/or when 
the coefficient c(x) varies appreciably. This is because a slight amount 
of predictor-corrector behavior is present in that an improved updated 
value of the pressure is used on the right hand sides of (76,77,78). 
(3) The computation time is much longer than the generalization to three 
dimensions of (68,69). 

In view of the experience (3) an implicit scheme is very desirable. 

3. Solution of the Diffusion Equation Implicit System 

In Section B an implicit scheme for solution of the diffusion 

equation was presented; the equations can be solved very efficiently 


as follows: Writing (57) in the form 





. n+l +B Dal qntl =D , m=1,....,M-1 
m mt+l m m m m-l m 
(80) 
it follows that 
oO oO 1 
= —_—o - = SS I — 
on 2( Ax)" Be (Ax)? a.’ = mn % 
(81) 
1 oO n oO n Go n 
= (—— - ee eee + == aco 
i re (Ax)? aa " 2(Ax)* eae 2 Ae Tal ° 


In the above relations dD, is assumed completely known at all positions 

m at time state n=o0. In addition, the conditions at x = o and 

x = MAx are either assumed known at all times or they are assumed known 
attime t=o and subsequently they satisfy a linear-discrete relationship 
equivalent to the usual Sturm-Louville problem. For the present, assuming 
the simplest case with homogeneous conditions (46) and (47), the system 
of equations (80) together with the left hand boundary condition (46) deter~ 


mine a one-parameter family of solutions, i.e., two sets of quantities 
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E and _ are desired which for any member ae satisfy the relation 
m 


oe eee ee (82) 
m m mtl m 


If (82) is true for any member it is true for m = o which requires by virtue 


of (47) that 


Eo = so, we = © £ (83) 
Oo Oo 

If now 
me = E a aes we (84) 
m-1l m-l om m-l 


is substituted into (80) the result may be written as 


Mert TT n a i or el 
ll eck oee ol 2 eos oe) 
m m m-l m m m-1l 
This in turn may be equated to (81) to obtain 
om 
SO = | as 
m m m-l 
mn 
Ernie | TE = CN yar la IE (87) 
m m m-l 


Starting with m=1 and the values (83), the En and Ee can be calcu- 

lated sequentially up to m = M-1. Once all the Ea! i are known (82) 
~~ 

allows the calculation of ae : sequentially in decreasing m, (m= M-1, 


+ 
M-2, ....-, 1) since form = M the temperature va? is known by (48). 


1 
This particular algorithm also has the advantage that the parameters 
A through en lie within a narrow dynamic range; it can be shown using 
the relations (81) that if En <1 then 
A A 


o SEBO SO c« = 1 
A 
m m MM 


Since by (83) EA =o , En lies between o and 1, Further, since the 
solution of the parabolic diffusion equation is analytic the solution we 
must be bounded and this with the bound on E requires that A also 
be reasonably bounded by (82). 

The greatest advantage of this particular implicit algorithm 
solution is the saving in computation; it requires only three multiplica- 
tions and two divisions per space step, per time step. It is a direct 
result of the matrix Ain (61) being tri-diagonal in form. The corresponding 
solution by direct matrix inversion, which would only need be done once 
if the coefficients are independent of n, would require in addition to the 
inversion procedure - the multiplication of a matrix by a column vector - 
which would be M multiplications per space step, per time step. 

4. Solution of the One-Dimensional Wave Equation Implicit System 

The above algorithm for the diffusion equation can be extended to 
the wave equation in certain cases, For the one-dimensional wave equa- 


tion (62), (63) a symmetric implicit scheme is 





_ q an 
m m es ie cel chal secant al 
At ne ee 
yntl " yn 
m-s m-s m Cc | n+l. n+] , 
At 2 Ax Pin Pm-1 Pal P ea g 
(89) 


To investigate the stability of such schemes the Von Neumann theory has 


been generalized by Kreiag- ante others, [11-15] 


Reference 4 summa- 
rizes the technique. Assuming that a suitable norm has been chosen for 


the domain of solutions of (88,,89) and that hence the Fourier analysis may 
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be applied to the system, let 


gt a 
m a k 

<i oma Hk) ek (m Ax) (90) 
Vv 5 k = —} © Vv 

m k 

a An+1 

m k 

_ jk(m Ax) 

i BS) an+1 | ° (91) 
Vv kK = = © Vv 

m k 


If (90) is substituted into (88, 89) the result may be written as (after 


cancellation of the common factor exp (j k( m A x) ) from all terms of the 











equations) 
A n+l oh 
Pk Py 
H ~ 4H = Oo (92) 
ar n+l 2 Y 
k k 
where Hy and A are the square matrices 
| cAt TOK 
il eae: sink 7 
H tk) 7 At SaCAt ok xX ea: Ax ares Ax 
Ax 2 s Does 2 
(93) 
mC Ake _. Ax 
1 : ae sin k 9 
H = 
>) Wii Herat aguls Ax sin ko aus Ax 
Ax 2 2° 2 
(94) 


The gain matrix G(k, Ax, At) is then defined as 
Gik, Ax, At) = (HH (95) 
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[4] 


and the general Von Neumann condition due to Richtmeyer may be 
stated as as follows: 
A sufficient condition for stability (and by the Lax equivalence 


matrices Hy and Ho is that the gain matrix 


o < Bt < 7% 
IG(k, Ax , At) | for o <nAt <T (96) 
all k in (90,91) 


be uniformly bounded. 


The matrix G for the present problem can be written as 





l-a°/4 ja 
lt 
G (kus < 2256) = ice a (97) 
1+a°/4 a l-a? /4 
where 
mee. 2c At : Ax 
ee k 0 (98) 
It satisfies 
G*G - GG*¥ = o (99) 
13] 


and is therefore unitary. Hence G is always bounded and the scheme 
(88,89) is unconditionally stable. The system may be solved by the same 
scheme as previously shown for the diffusion equation since it is also 
tri-diagonal. In order to guarantee the non-singularity of the matrices 
involved it is easier to deal with integer subscripts; (88,89) can be 


written as 


n n n+l n+l 


m-1 ‘m+tl — 1) (100) 
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n n n+l n+l 
= i : 101 
n+l on oy ee Ge) p Pt) ~ Pm-1? (101) 


Vv -v =- 
m m 4Ax 


The scalar equation (80) becomes the vector equation 


. untl +B yntl _c¢ yntl me. (102) 
m —m+l i an m—m-1 — 
where 
n+] nD erat oe Laan 
yntl Pin Pm 4Ax m+1 m-1 
= = , dm = (103) 
00 n+l = yt : c At ce a 
m m 4 Ax m+l m-1l 
and 
cAt 1 . 
o 4 Ax 
= = , =-A , 104 
c At “ : 5 1 : ue) 
4 Ax 
Again letting 
u = E u + f mn — Ogle. ses, M (105) 
snes nee mM 
the results 
ee Ce ee ee ae OG) 
m m m m-l m 
-] 
a 7 Bo : Cine in Ctmer) Male weer Mol 
(107) 


are obtained. To solve the system Suppose the matched boundary condition 


is applied at x = MAx inthe form 


p eel p 
m-1 m 
= (108) 
V ih O V ® 
m-1 m 
and a homogeneous velocity condition at x = o in the form 
p p_(n) 
° |) = = , P (n) given. (109) 
O 
vs O 
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The "matched" condition (108) is clearly only an approximation to the 
correct condition for large Ax, but becomes progressively better as Ax- 0. 
This closed=trough problem has been programmed on the IBM 360 with 

the following results: 


A, The computation time for equivalent Ax is reduced by an 
order of magnitude over the scheme (73,74) 


B. The accuracy for a given Axis considerably poorer than 
obtained by either (68,69) or (73,74) but this is thought to be 
because the boundary condition (108) is properly simulated 
only in the limit of Ax - 0. Because of (A) above a much 
more rapid refinement of Ax is possible leading to very 
satisfactory results. 


Variable parameters are routinely handled and no instabilities 
were observed for reasonable parameter variations. 


Q) 


59. Three-Dimensional Implicit Formulation; Fractional Step Algorithm, 


The results obtained for the one-dimensional implicit scheme for 
the wave equation, particularly the increase in computing speed, make a 
three-dimensional implicit scheme highly desirable. Therefore, a 
straightforward generalization of (100,101) was considered. The resulting 
scheme is 
p(2,m,n,rt¢l) - p (£, m,n,r) 


Ee - Bist mana) ts (v ( £+1,m,n,r) - (£-1,m,n,r) 


tv (441 ,m,n,rtl) - v_(2-1,m,n,rtl) 





as v (£.m4] ,n,7) a Bee ‘ m-1,n,r) af vit ‘ mae) , ntl) 


1 
Ay 
=v ,m-l,n,rtl) ) 


+ (v (4, m,nt1,r) - v (4 ,m,n-1,r) +v (4 ,m,ntl,r+l) 


Az 
-v (£,m,n-1,r41)) | (110) 
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ae anginynt1) - v4 anyanyn) i 


ee IS n) 


4 Ax (o( 2 a mynyb), Seek L.~1 gmenyr) + p(b+1 man, tt!) 
~p(£-1,m,n,rt+1) ) (111) 
v (£,m,n,r+1) -v (2£,m,n,r) = 
y! vs 


a (p(2,m+l,n,r) - p(2£,m-1,n,r) + p(Z ,m+1,n,r+1) 


- p(£,m-1,n,r+1) (112) 


v (2 ym,n,rtl) - vi (4 ,m,n,r) = 


Atc(£,m,n) 
- “4 =  (p(£,m,n+l,r) - p(£,m,n-1,r) + p(£,m,n+l, r+1) 


~ p(2,m,n-1,r+l) ) (113) 
where as before p(£,m,n,r) = p(LAx, mAy,nAz,r At), etc. This 
algorithm is again unconditionally stable as G(k, Ax, Ay, Az, At) is 


unitary but the tri-diagonal behavior has been lost; the matrix H, is of 





il 
the form 
sin k_Ay, j c Sinek Az | 
TOA a= Me 3 
: ] , O : O } 
Hy At 
Oo 4 1 F O 
O ’ O ’ ] (114) 





which is penta-diagonal. Several problems were solved using this algorithm; 
it is highly accurate, stable under quite general conditions, but is relatively 
slow and requires a great deal of programming effort to use in the applica- 
tion being studied. 

Many investigators have proposed schemes for solving the three- 


dimensional wave equation efficiently while maintaining stability. The 
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Li 4) (17) 


Soviet mathematicians Bagunovskii and Goduno, and later Yanenko 
have used partitions of the spacial operator for the various derivatives, 
applying each operator for fraction of the total time step At. Douglas and 


[18 ] 


Gunn proposed a very general scheme for any ‘number of space 
dimensions which involves p-tridiagonal calculations for a p-th order 
cartesion system. Douglas and Gunn's paper contains examples for the 
diffusion equation; the procedure has been applied to the system (9,10,11, 
12) of Chapter II under the symmetric condition Ax= Ay= Az. The resulting 
system of equations -- which is quite long -- is shown in Appendix A, 
The system can be shown to be unconditionally stable if the coefficient c 
is a constant; for variable c considerable computational experice seems 
to indicate considerable margin of stability exists if the variations in c 
do not exceed five percent, 

6. The Fast Fourier Transform Technique for Constant Coefficients 

In the special case that the velocity c may be treated as a 

constant over the domain of integration, the use of the digital discrete 


9] 


Fast Fourier Transform (FFT) algorithm provides a simple, elegant 
solution of either explicit or implicit algorithm schemes. The application 
of the discrete Fourier transform to this problem in the manner to follow is 
new. 

Since implicit algorithms in general, as shown above, exhibit 
Superior stability properties and when properly formulated greatly reduced 
computation requirements, the FFT algorithm was applied to their solution. 

The implicit scheme (100,101) will be solved to illustrate the 
method. Assume that the initial states 
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V 
m 
are known and that the number of sample points M, including boundary 
points is an integral multiple of two. This restriction is not absolutely 
necessary but its imposition guarantees the maximum computational 
efficiency for the FFT algorithm. The discrete set of values may then be 


represented by the finite inverse Fourier transform as 


xO 
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Pin ae Px en k 
= eed e! M ” an) 

O k=o AO 
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A 
where the spectral functions Py a are obtained by the process 

Ko oO 
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This result, coupled with the fact that, by (116) 
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, = > i e M , P= integer, 

ea K Ome (118) 


leads to the following equation if the transform is applied to (100,101) 


under the assumption that c is a constant, 
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This can be true for arbitrary m only if the large bracket term vanishes, 


Hence defining the matrices 
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G, (k, AxXpodt)e = (120) 
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1 Oo 
G ({k, Ax, At) = (121) 
O 
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PWHESt 21k 
a= jonx I9 (122) 
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Ph Mol Py jee m 
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Hence apart from the calculation of the powers of G 19 *) the 
entire process involves only the application of the FFT algorithm at the 
time states n=o and the desired state n . The powers of Gy (k) may 


be most efficiently evaluated for large n by use of Sylvester's 


[20] 


theorem, 


The procedure is readily extendable to three dimensions; the 
penta-diagonal scheme (110,111,112,113) has the solution 


pi( Lean, npr) 
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There is, of course, an unfortunate oversight in the above manipu- 
lations: the equations (100,101) involve the boundary conditions at m=o 
and m=M_ which are outside of the indexing of (100,101). The same must 
be true in the above formulation so that what has been solved is the system 
in which p(£,m,n,0,),v(2,m,n,0) are known and the tacit homogeneous 
boundary conditions p(-1,m,n,r) = p(L,m,n,r) = p(£,-1,n,r) = p(1,M,n,r) 
=o, etc., have been assumed, (A homogeneous nearest-neighbor shell 
of boundary conditions on p,v enclosing the parallelepiped LAx by MAy 
by NAz.) The method of solution can be extended to handle this case by 


[31] 


making use of the theorem by Courant and Hilbert, wherein it is shown 


that: homogeneous differential equations with non-homogeneous boundary 
conditions are essentially equivalent to non-homogeneous differential 
equations with homogeneous boundary conditions, Usually the applica- 


tion of the theorem requires in addition that certain smoothness conditions 
be satisfied by the boundary conditions, such as differentiability, etc. 
The theorem is a natural extension to partial-differential equations of the 
familiar result that the solution to an inhomogeneous differential equation 
consists of a solution to the homogeneous equation which satisfies the 
same initial conditions plus any particular integral of the inhomogeneous 
equation satisfying homogeneous initial conditions. The theorem is 
restricted to linear equations. 

To apply this theorem to the problem under discussion several 
methods are possible: (1) the application of the Lagrangian multipliers 
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method, (2) the eigenfunction expansion method in which the projection 
of the boundary conditions on the complete set of eigenfunctions 
appropriate to the linear homogeneous operator equivalent to (100,101) 
is used to augment the operator, (3) the use of symbolic functions to 
include the boundary conditions. 

Method (3) is the easiest method to apply in the present case. 
The preliminary details of the method proceed by finding the adjoint 
operator A* for the system assuming homogeneous boundary conditions; 
the scalar product 

(A*u, A, = (A*, A, u) (129) 

is formed where A,is an extended linear operator and the domain of u 
now includes the non-homogeneous boundary values. The forced equality 
in (129) generally requires that A, involve singularity functions; in the 
continuous differential equation case usually the Dirac function ei its 
derivatives are involved. In the case of discrete systems the Kronécker 
function 


oe tl aa 


= Oo eas (130) 
replaces the Dirac function, 


For the problem under consideration here let the solution at time 


state n=o be denoted as ue : 
wu = uw (mAx) 
O 
p (m Ax) (131) 


vo (m Ax) / . 


O . O ; 
u_is assumed to consist of two parts: a:homogeneous part u which is 
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the solution to the linear difference equations (100,101) with 


u (m) ' mo or M-l 


= Oo : m=o or M-l, (13:2) 


u, (m Ax) 


s O & e e 
and an inhomogeneous part uu. which satisfies 


u? (m Ax) = 0 , Mm? oF Or ll 
=i 
=u (0) ae oe, 
i 
O 
= u. (M-l) , m= M-l ., (133) 


The values ue (o), ue (M-1) are the boundary values for the problem and 
are assumed knowna priori. In terms of equations (100,101) the predicted 


value ti has component ! - satisfyin 
me p Ba se ying 


1 .o _ cAt ae uel Pee 8 \+6 a: ae 
Pin Pn 4Ax m+1 m= m+l m-l m,oPo m,/M-1?M-1 
(134) 
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(135) 
For the applications which require a linear combination boundary condition 


such as for example 


O oo 
Po OY ait] (136) 


where b is a constant- the matched boundary condition is an example - 
O O 
g replaces p- in (134). 
Now denoting the FFT of ue as ne and of fa as He and applyin 


the transform to (134,135) results in the matrix equation. 


_-j_27 (M-)) k 


A 
Gu. = G u + uo) +e M uy (M-1).. (137) 
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By hypothesis uy is zero for all m in (134,135) as far as the terms 


O 


involved in G, u, are concerned so that 
SO 
Gu.=o (138) 
Cy mea 
and hence 
_ = Gi { u” (0) + Pe i 5 w(m-1) } 
oil i: = =a (139) 


Al 
is the solution for the FFT spectrum of the inhomogeneous solution u, at 


~ 


] 
time state n=l At. The solution for u 


h is that already found in (122) 


1 
h 1 or ia 


Ic > 


A 
= Gig u : (140) 


Since the FFt is a linear mapping operation the spectrum of the 


sum of two functions is the sum of the individual function spectra and 


hence 
Al | A] 
u = u + u. 
= =i 
AO -l] ao 
+ 
Gio UL G, ui (141) 


The solution (141) represents the spectrum at time state n=l ,. To 
proceed to time state n=2 the process is again repeated except the homo- 
geneous spectrum at time state n=l is now taken to be the total solution 


spectrum given in (141). The equation (141) for time state n=2 thus 


becomes 
f= ogi! + of! a 
ill, + a) 8) +) by 
“Gy Prager ees. 42) 


where 
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al -1 1 Kk i } 
— ~ = 
| 124 ene " ee (143) 
By induction the general solution at time state mAt is given by 
‘Ts nao n-l_-lao 2. -lal -lan-1l 
u =%G,,.) u, + (G, 9) G, u, +(G,,)  G, u;, t..0.G ui. 
(144) 


An obvious thought is that since the matrix G is unitary and bounded 


10 
and hence the homogeneous solution is stable independent of At, why 
not set At = t the desired final time and iterate the process only once? 
The answer to this conjecture is negative. The system of equations 
(134,135) has a stability criterion which no longer is unconditional, 
depending in a complex way on the functional behavior of the boundary 
values. A moments thought leads to the realization that this must be so 
Since only the homogeneous solution is effectively implicit; the inhomo- 
geneous solution by (139) is explicit. Hence, the overall scheme is only 
semi-implicit and the magnitude of the At allowed depends on the magnitude 
of the increments in the boundary conditions between time steps, 

Since the rate of change of the boundary values with time is 
linearly related to the frequency contained in the frequency spectrum of the 
boundary values, and as will be shown, the magnitude of the frequency 
must be low to guarantee proper non-dispersive simulation relative to the 
Shannon rate, the increments in the boundary conditions are generally quite 
small. This, coupled with the fact that only the inhomogeneous portion of 
the solution is not unconditionally stable, leads to the qualitative specu- 


lation that significant magnitude At may be used. 
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Computer solutions of the above algorithm have verified that such 
is indeed the case; quite stable and accurate results were obtained for 
At's more than twice the size of that allowed by the scheme (71,72) and 
with reductions factors in computation time of greater than ten. Since 
the number of calculations for the three-dimensional problem is of the 
order of the cube of that for the one-dimensional problem, the use of 
schemes like (110,111,112,113) is only justified for problems in which 
the variation in the coefficients is so rapid that prohibitively small domains 
would be necessary in order to use the FFT scheme. 

An additional advantage of the FFT scheme is that for problems 
of the type of primary interest in this work in which spacial partitioning 
is possible, the objective of integrating the wave equation is primarily 
to find the steady state generated values of the unknown boundaries of a 
parallelepiped in terms of known values over the other boundaries of the 
parallelepiped. (For example if the values were known on six surfaces of 
a parallelepiped, what are the steady state boundary values on the seventh 
surface?) By virtue of (126), this would correspond to inverting the 
transform only for a fixed two-dimensional plane so that only L M-point 
inverse transforms need be calculated rather than (LxN) M-point inverse 
transforms - a considerable saving of computational effort. 


D., BAND-LIMITED PROPERTIES OF DISCRETE MODELS;: SPACIAL SAMPLING 
DIMENSIONS 


The above algorithms, whether continuous-discrete, explicit-discrete, 
or implicit-discrete, all require great amounts of numerical computation, 
the amount depending fundamentally on the spacial partition parameter 
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magnitudes Ax, Ay, Az. In section Ait was claimed that the element 
partition size Ax, in the one-dimensional case, was determined by the 
operating frequency whose simulation is to be studied and the sound speed 


C. 


i. CutotimeProperties 


The simple discretized wave equation (68,69) is an exact dual 


of the one-dimensional periodic, (in the constant coefficient case) mass- 


5] [23] 


loaded line originally studied by Baden-Powell, Lord Kelvin 


others. The propagation properties of such lines can be found by appli- 


[24] 


cation of Floquet's theorem; if the assumed solutions 


p(mAx,t) = Ae eee (145) 


v(mAx,t) = view" sede ese) 


are substituted into (68,69) the determinental equation for the propagation 


constant syaresults, 
2 
WwW 
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( sinh kop = - (oamayie . (147) 


The form of the result (147) varies slightly depending on the finite difference 
formulation used but all discrete schemes result in the same Qualitative 
behavior. For lossless propagation y must be pure imaginary 

y= 5B, (148) 


which gives in (147) 


sin PAx = + = Ae 
. : oe (149) 


Real B solutions to (149) occur only if 
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a za (150) 


from which the cutoff frequency be is obtained as 


Cc 
ea 151 
2 27 Ax (151) 





Since the velocity c is an inherent property of the medium, if propagation 
of a monochromatic signal of frequency f is to be simulated, f should 


satisfy the inequality 


i < 2 oe : when 


Assuming a nominal velocity in sea water of 1500 meters/sec. this 


restricts x to be less than 


Ax < _ : CESS) 


For a frequency of 500 Hz., Ax < .478 meters, a very short distance. 
The Shannon spacial sampling distance would be 1.5 meters so that the 
cutoff frequency restriction is more severe than the sampling restriction 
necessary to properly simulate the partial-differential equation. 
2. Dispersion Properties 
Equation (149) also indicates that the propagation is dispersive, 


i.e., the propagation constant 8 , 
(154) 


is not a linear function of frequency except for very low frequencies. In 
order to make the argument of the arcsine in (154) small, Ax must be 


restricted to 


- el (155) 
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or 
pif (156) 


for non-dispersive simulation. Thus the restrictions on spacial increments 


are very severe, especially for very higher frequencies. 
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IV, COMPUTER RESULTS FOR THE MATCHED BOUNDARY CONDITION 


A. PARTITIONING OF THE SPACIAL DOMAIN 

The integration of the acoustic wave equation over large spacial 
domains in which the sound speed c varies with position is usually 
impractical to carry out on digital computers because of the hug amounts 
of core required for even crude finite-cell discretizations. To implement 
for example the algorithm (110,111,112,113) for a rectangular parallele- 
piped region consisting of (25 x 25 x 10) spacial integration volumes 
requires a minimum of 50,000 words of memory just to store the required 
values of the variables. The coding of such a problem requires great 
amounts of indexing calculations which typically would require at least 
an additional 10,000 words of memory for the stored instructions , and 
temporary work storage. These requirements, especially in view of the 
fact that the parallelepiped would simulate a quite small physical volume 
for all but very low frequencies, has made ray theory and normal-mode 
model approximations the only practical methods being Ease employed 
for investigating undersea acoustic propagation. 

A problem of considerable interest is the calculation of the steady- 
state acoustic loss over great distance for various frequency acoustic 
waves in seawater wherein the sound speed c varies with position. As 
pointed out in the introduction, at present this calculation is carried out 
using ray theory to establish an approximate sound trajectory and then the 


loss is determined by integrating an empirical loss factor per unit length 
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along the trajectory. 

In view of the practical limitations on core total storage size for 
available digital computers, the more fundamental direct integration of 
the wave equation is practical only if problems involving large spacial 
domains can be partitioned into relatively small spacial domains which 
Can occupy core storage sequentially. The macro domain of 25 x 25 x 10 
micro domains, (micro domains will refer to fundamental small elements 
of size Ax by Ay by Az) might constitute such a relatively small 


spacial domain. 


B. THE INHOMOGENEOUS LINE IMPEDANCE 

In one-dimensional uniform transmission line theory a finite section 
of line of length £ can simulate a portion of an infinite line provided 
the voltage and current at the ends of the line are constrained by the 
equations (see Fig. 3) 


v(o,t) = Z. i (o,t) (157) 
yale Cail liao Zi L,t) (158) 


where a is called the characteristic impedance of the transmission line 
and is related to the series impedance per unit length Z and shunt 


admittance per unit length Y, by the relation 


h 


S = : (158) 


If the line is lossless and possesses only series L and shunt C 
elements the equations of the line are 


C ae ee (159) 
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Oo. ene 160 
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which are the complete duals of the equations (62) and (63). If steady 
state solutions are obtained for (159,160) by assuming 
A 
v (x,t) V (x) . 
ae a (161) 


i(x,t) (x) 


A 
the voltage and current phasors V(x), { (x) satisfy the relation 


a ee (162) 


where 2(x) is the complex impedance 


Z(2) cos(w/LC(£-x)) + jZ. sin(w/LC (£-x)) (163) 
O Z. cos (w/LG (£-x)) + j2() sin(w/LG(£-x)) 





ZA(x) = 2 


and 
a = ae 
oh 
L 
= C 3 (164) 


Thus the characteristic impedance is a pure real constant only if the 


line is lossless or in the happy circumstance where 


J ace (165) 


in which case though a and Le are both complex their phase angles 


h 
are equal and thus cancel. In (163), 2(9 is the terminating impedance 
of the line; if 2( has the value 3, (163) gives the result 

Z(x) = ae ; (166) 


a pure constant. The general steady state solution to (159) may be written 


as 
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¥ (x) = Vie ssh + Ve ie sh : (167) 


if the line is lossless this results in 


Wave ON ye ee (168) 
f (x) = Le (169) 


where ve ; V are the amplitudes of positive and negatively traveling 
waves respectively. Causally if a monochromatic disturbance is propa- 
gated in the positive x direction with amplitude ve the coefficient 

V will be zero unless there is a change in the characteristic parameters 
of the transmission line. In such a case reflection coefficient R can 


be defined by 
yO Jwlls -x) 


ye qi @vLC (2 -x) 


R(L) _ jew J LC (£ -x) 


R(x) 


(170) 


a 


and is a measure of the "non-match" of the line, If the parameters L and 
— (25a 

C are not constants but functions of position Shelkunoff has general- 

ized the concept of reflection coefficient for inhomogeneous transmission 


lines in terms of the propagation constant y (x) where, 


Z (x) Yan (x) 


II 


¥y (x) 


a (x) +jB(x) . (171) 

The derivation to follow deviates slightly from Schelkunoff's approach 
being specifically oriented to the acoustic problem at hand. Consider the 
two semi-infinite half spaces separated by the plane x=o, The region 


x < o is assumed to be homogeneous with constant parameters c, p, x, 
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the region x > ohas slowly varying parameters x(x) and c(x). Assuming 
steady state dependence of the form 
A 
p(x,t) p(x) 
ae | , e) ue (172) 
v (x,t) v(x) 
in equations (62,63) reduces the partial differential equations to the 


ordinary differential equations 


ar a = 1 . 

= jwpv (173) 
dv 

—se — ee —-_ 174 
ax ee rd ( ) 


For x <o the solution to this system for P(x) assuming unit incident 


amplitude is 


B(x) = e **4R(x) e ” * (175) 


where R(x) is the reflection coefficient. Since the reflection coefficient 
at any point Xo is a measure of the mismatch as it affects the behavior 
for x < Xo of inhomogenieties occurring for x > Xo (175) as pointed 
out by Schelkunoff, provides a plane wave reflection coefficient which 
holds for any x, The propagation constant y is obtained from (173) 


and (174) as 


7 P 
_ ,_W 
j Gc) : (176) 
Defining the general impedance at a point as 
A 
Bx) = —ed (177) 
V (x) 


and using (173), the result may be written as 
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4(x) le 
~~ Pyar 
jwo \dx 


= LWP (178) 


From (175) 


B(x) =" lan 


-y (e ’*_R(x)e” *) + R(x) e” * 


-y X yx 
=~ J Wo e + R(x) e (179) 


is eet R(x)e” * 


if R(x) < < 1 and slowly varying. By virtue of the continuity of the 


A UA 
ratio P/V the function (179) is continuous across an interface boundary 


even if x(x) and c(x) are discontinuous. By virtue of the fact that 
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A d j A Z A a? i A pn 
= Tae (log P) ) + P —y2- (log P) -y P (182) 


or 


Zz 
d A 
Py (log P) = y 


dx 


*- (og?) * (183) 


By use of (178), equation (183) becomes 


G_ ,; =jwo, _ ,4 ,;r-jwp 

dx ( Z (x) y | Z(x) ), oe) 
or 

dz yilce 

dx" jwp PO _ 


This is a first-order non-linear equation for the impedance as a function 
of the propagation constant y (x). Alternately, (179) may be inverted 
to R(x) in terms of 2(x) and the result substituted in (183) to obtain 
the differential equation for the reflection coefficient 


dR 


eS —(1-R”) = loonie) B= 208i = On (1 86) 


a first-order non-linear equation of the Riccati type. If 


S Gy fog GW) << i (187) 


then (186) becomes 


dR 
— Lit gs 2yR (188) 


which has the solution 


Xx 
-2 
R(x) = RO) e Jy (§) dé, : (189) 
By virtue of (176) it can be seen that the effect in this case is simply a 


phase shift due to the length of path and to first order if c(x) is smoothly 
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varying, (in a multi-dimensional problem) to cause a bending in the 
direction of the wave. uj 
C. THE MATCHED BOUNDARY-CONDITION DOMAIN 

The matched termination or characteristic impedance has been exten- 
sively discussed in the excellent book by Brillouin 5]. he considers the 
general periodic structure lattices consisting of polyatomic molecules 
and derives an expression giving the characteristic impedance as a 
function of the periodic cell internal energy. He concludes his discussion 
by emphasizing the difficulty involved in theoretically defining the 
characteristic impedance for structures in which more than one dimension 
is involved or where more than nearest neighbor interactions are 
considered. 

Here the objective was to partition variable-parameter, large domains 
in three-dimensions with the objective of computing steady state 
trajectories and pressure distributions. Consider Fig. 4 which shows a 
large domain consisting of 27 macro-domains. Each of the individual 
macro-domains consists of micro-domains; (L=M=5, N=3). As the 
simplest and most straight-forward method to apply, an Adams-Moulton 
predictor-corrector, with dynamic error control program called RKAM (see 
program listings), was used to integrate the explicit three-dimensional 
finite difference scheme equivalent to the one-dimensional scheme shown 
in (68,69). The procedure for integrating the large domain was to integrate 
the problem two ways and then to compare the results: 

(A) Integrate all (3+1)(5)(5)(3)(27) = 8100 equations in an overlay 


program until total steady state energy was attained for the 
entire large domain with characteristic boundary conditions 
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being employed on the outer-most boundaries. The energy 
source was placed asymmetrically in the interior of the 
innermost macro-domain in order to preclude any possible 
cancellation effects due to symmetry. 


(B) Integrate the 300 equations for the innermost macro-domain 
with the energy source located in the same position as (1) 
above. Characteristic boundary conditions were again 
employed except their location was the boundary of the inner- 
most macro-domain, Upon reaching steady state conditions 
the FFT was employed to find equivalent oscillators for the 
pressure and three velocity components at each micro- 
domain on the surface of the inner-most macro-domain. These 
oscillators are then assumed to be the sources for the nearest 
neighbor adjacent micro-domains in the adjoining macro- 
domain, (the boundary micro-domains are considered to belong 
to both adjoining macro-domains.) In general, only the 
nearest neighbor macro-domains are considered in establishing 
steady state values for the boundaries. Since the medium is 
linear, the rms magnitude of the states on the boundaries can 
be found as the linear superposition of the quadratic contri- 
butions to the rms value of each adjacent macro-domain 
boundary effect. Proceeding macro-domain by macro-domain, 
the equivalent rms boundary values for the large domain were 
obtained. 


The above procedures were performed on the Fleet Numerical Weather 
Facility CDC 6500 computer which carries out floating-point calcula- 
tions to 14 decimal digits. In the case of constant coefficient domains the 
above two procedures produced results in excellent agreement considering 
the prodigous amounts of computation involved; the two computations were 
in agreement at all but a very small set of values to three decimal digits, 
There were no major discrepancies, 

In the variable-parameter domain case the concept of matched charac- 
teristic impedance looses significance unless more than nearest neighbor 
interactions are considered because the characteristic impedance can be 
fundamentally defined only for periodic domains. However, when the 


variations in the velocity did not vary over one percent between adjoining 
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macro-domains the agreement of the outer boundary values was still better 
than one percent. 

An approximate procedure to take into account the variations in c 
immediately beyond the boundaries of a macro-domain was investigated, 
This involved using the solution to the differential equation (185) to 
improve the value of z (x) to be applied at the boundary micro-domain. 
Since only the nearest neighbor micro-domains should effect the local 
impedance on the boundary, the impedance at the micro-domain two 
micro-domains into the adjacent macro-domain normal to the boundary 


was assumed to have the value 


a(2 dx) = ,| o% ? Ax) (190) 


Equation (185) was then numerically integrated back to the boundary to 
obtain Ao), (the improved value to use for the ratio of P (bnd)/V (ond) ). 
With this added refinement the boundary values again were in agreement 
to three significant figures for variations in c of as much as five per- 
cent between adjacent macro-domains. Variations in c as large as five 
percent seldom occur over macro-domains as small as considered, 
1. Stability Of The Internal Energy Criterion 
The criterion used for steady state determination in a macro- 

domain was the total amount of internal energy contained in the macro- 
domain. The kinetic energy of a micro-domain is given by (at any time 


state r), 
= + pax Ay Az lv (4,m yn x) + v (#,m,n,r) + v) (£,m,n,r) | 
micro 2 ae ee yaa Zi oil 


(191) 
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and the potential energy is 


ee ee (192) 


V (2,m,n,r) = rn 


iilcre 


Hence the macro-domain has the energy 


= aT + V 
oes u ( micro micro. (193) 
i, Vial) 
For very low frequencies the function U varies instan- 
macro 


taneously but periodically with the time at twice the frequency of lowest 
harmonic present in the medium. Thus the criterion for steady state used 
was to monitor the instantaneous stored energy every period of the lowest 
frequency being simulated. This value was found to be stable to five 
decimal places. 
2. Transient Behavior Of the Models 

The behavior of the various integration schemes to a step input 
is illustrated by the simple model (72,73). In Appendix B is shown the 
step input pressure spacial distribution response of a 100 element one- 
dimensional domain of unit length ( Ax=0.01). The responses are 
shown for t=0.,0.2,0.4,0.6,0.8,1.0,1.2,1.4, and 6 seconds. The 
speed c was taken to be unity. The symmetric form of the equations, 
(20,21) were programmed so that the characteristic impedance was unity. 
The results show the impulse excited transient oscillations characteristics 
of such models; however the matched ccndition is seen to produce no 
distinguishable reflections and predicied velocity of unit/sec is quite 
accurately simulated, The internal energy reached steady state to five 


decimal places in three seconds. 
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In Appendix C is shown the result of applying a unit amplitude 
cosine wave of frequency 1 Hz to the same line as in Appendix B. Again 
the characteristic transient oscillations occur along with some minor 
waveform distortion in the region of the maxima and minima of the cosine 
wave, These oScillations are caused by the shock excitation of the 
micro-domain self resonant frequency. The waveform distortions persist 
for a considerable period of time even after the energy has stabilized on 
the model, This emphasizes that caution is necessary in deciding when a 
model has reached steady state. Inthe problem shown an FFT of the 
spacial distribution could be used to indicate that only the desired 
simulation frequency is present but for three dimensional problems sucha 
procedure is probably computationally too extravagant. 

3. Two Dimensional Model 

An example of a two-dimensional model integration with variable 
velocity c as a faction of one of the dimensions is shown in Appendix D. 
The program used was WV2525, (see program listing). An initial condi- 
tion of pressure equal to unity for x=o and y varying from o to 25 Ay was 
applied. The velocity for y=25 Ay was three times the velocity at y=o. 
For stability the integration step size required was very small, Amite 
sec, The two-dimensional continuous-discrete system equivalent to 
(68,69) was used, 

In Appendix E the small amplitude transient response of the 
plucked membrane problem is shown, This is an example of a non-zero 


initial condition natural response simulation. 
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V. APPROXIMATE FAST FOURIER TRANSFORM COEFFICIENTS 


A, INTRODUCTION AND DEFINITION OF THE TRANSFORM 
The Fast-Fourier-Transform (FFT) algorithm provides an efficient 
digital computer algorithm for performing the one-to-one discrete 


sample space mapping, 


_ 27 (n At) 
| aoa Te 
G(k) = N ~ g(ndAt) e (194) 
n=o 
N-1l 2a K 
WAT t 

g(nAt)==Z Gk) e Be, ‘ (195) 

k=o 


where g(n At), G(k) are the object and image sample spaces and At is the 


sampling interval, The properties of these transform pairs have been 


L19,26,27] For monochromatic 


extensively discussed in the literature, 
signal sequences g(nAt) the spectral function G(k) can be guessed or 
arrived at by heuristic reasoning but for general functions the summation . 
in (194) can usually only be carried out analytically provided the parti- 


cular finite series have known sum forms, 


The continuous function Fourier series coefficients a where 


co eat ; 
g(t) = x ome (196) 


=- @ 


can be obtained either by integration which is the continuous process 


[6 J 


analogous to the sum (194) or by means of symbolic differentiation 


8] The possibility of using the symbolic 


using impulse functions, 
method - modified for use on discrete sample spaces - is the chief result 


of this chapter, 
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B, REQUIREMENTS TO OBTAIN COEFFICIENTS BY THE IMPULSE METHOD 
In the case of continuous functions the requirements for applicability 
of the impulse method of finding Fourier series coefficients are: 
(A) A spectral representation for the Dirac function 6 (t). 


(B) An analytic expression relating the derivative Fourier 
coefficients to the original function coefficients, 


(C) A piece-wise analytic expression for the function to be 
represented. 


(D) Atheorem giving the derivative for discontinuous functions, 
It will now be shown that when requirements analogous to the above are 
satisfied for discrete functions, the FFT coefficients G(k) can be 


obtained, 


C. SPECTRAL REPRESENTATION FOR THE SYMBOLIC IMPULSE FUNCTION 
The symbolic impulse function for discrete sample spaces is the 

Kronecker delta; it is defined here as the mapping function 6 ny 874) 

-1 


B (5-4 f(s) = £9 (197) 
O 
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where f(s) is a suitable testing function defined on the discrete set of 
points s=o,...., N-l andis periodic modulo N; i.e., 

f(s) =f(stN) , alls . (198) 
This definition is equivalent to the usual definition of the Kronecker delta 
as 


6 (s-2) = ] ; Sinagall} 


=o : Se ae (199) 


The function possesses the FFT coefficients 
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_onk 


ed N ‘ (200) 


Glk) = ee 


as can be seen by direct application of (197) and (194). Hence 


N-l , 27k (s-0 
6 (s-9 yee : (201) 
k=o 


D. DEFINITION OF THE DISCRETE DERIVATIVE AND ITS SPECTRAL 
REPRESENTATION 


Many choices are possible for defining the discrete derivative of a 


discrete function; the one used here is the backwards derivative commonly 


defined in the literature, 9] 
, O(n At) preiiin=1) at) 
g‘(ndt) = SA! me arta (202) 


It is a linear definition possesing the inverse 
g(nAt) = g((n-1) At) + g’ (nAt) dt , (203) 
or predictor 
Nae 
g(nt At) =gein- Aare gq (s AtpAt, (204) 
s=n 


where o<n Ss N-l. Using this definition in conjunction with (195) 


results in the derivative spectral representation 











27k 
-j 2 
a tine UN j Ta ae 
G nat) =e Gh) == _e (205) 
k=o a 
or in general for the mth derivative 
a Fa OTK 
N-1 _ N me j (n At) 
gmat =5 Gk) flier } Sear (206) 
k=o 
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An important result is obtained by applying (206) to (200) 


oe m .27k 
= { l-e = IN At (s-9 At 
— } e (207) 





yt? (s-2)'At)& 
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E., ANALYTIC ENVELOPE EXPRESSIONS FOR DISCRETE FUNCTIONS 
An analytic-envelope representation for g(nAt) is necessary in order 
to be able to carry out discrete derivatives according to definition (202). 


Typical examples of such representation are 


(A) g, (nAt) nAt , osn S N-1 


(8) g, (nt) inAtl |. «UG hn em 


lA 


(N/2) - 1 


n At , o=@ 


(C) g 3 (n At) 


= -nAt. , N/2 <n < N-I1 


27 
N At 





(D) g,(n At) = sin ( (nAt)),o <N-l. 


The functions (A) and (B) are polynomial type functions with no "discontinu- 
ities in analytic-evelope" over the entire interval o <n s N-1. (C) is 
an example of a function described by two different envelope representations 
within the intervalo <n s N-1; it is said to have a discontinuity in 
analytic envelope in the region d n= N/2,. (D) is an example of a tran- 
cendental function representation. 
F, THEOREM ON DISCRETE DERIVATIVES OF DISCONTINUOUS ANALYTIC 
_ ENVELOPE FUNCTIONS 
Aaaloaonet to continuous function symbolic theory, the unit step 


function Usy (s-£), of a periodic sequence of discrete sample points s=o, 


1, «eee, N-l is defined as the linear functional mapping 


78 


-] N-1l 
Us) (s-9 f(s)= Ss f(s) .« 
O Ss=L 


(208) 


9, 44 


From (194,195) itis clear that the results for the coefficients are indepen- 


dent of At; it thus will be dispensed with in subsequent analyses, 


Us, (s-2 as defined in (208) is equivalent to the piece-wise definition 


Us, (s-4 =] , Cae 2 
=o 2 Gi Le (209) 
This function possesses the property that 
N-1 N-1l 
a Us. (s-9 ny i= a (Us (s-9 ~ Us, (s-1-9 } £(s) 
s=£ S=O 
N-1 N-1 
es ral 
=f s=L +1 
= f(£ ) 
N-l 
ners eel), (210) 
S=O 
so that symbolically 
(211) 


Us, (s-9 = 6 (s-9 
A change in a continuous analytic-envelope function can thus be described 


by 


g(n) = g, (n) +O Us, (n-2) (212) 


where g, (n) has the same analytic-envelope representation for all n. This 
represents the situation of a discontinuity in analytic -envelope at the 


point n=. See Fig. 5. The discontinuity is a jump of magnitude a. 


If (211) is used in (212) the result for the discrete derivative is 
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g‘(n) = g, (n) ta6 (n-£) . (213) 


This proves the theorem: 


Theorem: The discrete derivative of a discrete function possessing a discon- 
tinuity in analytic-envelope at the point n=£, is the backwards derivative 


defined in (202) for all n #2 plus a Kronecker delta function at n=2 of 


amplitude _q@_, where 


a= g(%) - g, (4-1) . (214) 


G. EXAM PLES 
1. Coefficients of a Sawtooth Waveform 
The preceding results and theorem can now be applied to a simple 
problem to illustrate the method. Consider the periodic sawtooth wave- 
form shown in Fig. 6 (a). The function has the representation 
g(n)igeeen , on <N-1 
g(n +N) = g(n) . (215) 


Assuming g(n) has an FFT then 


_ 27k 


| N-1l j 
ain) = eee (eee 
k= 


O 


i 


=n, - o <n s N-l, (216) 
If (210) is formally differentiated by means of (196) there results, (see 
Fig. 6 (b) ) 


n- (n-1) 


g ‘(n) 


= ] ; l sehen (217) 
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and 


g * (n) o - (N-1) 


=e ; n=o . (218) 
The results (211) and (212) can be combined into the single result 


g’(n) = 1-N6,) , osn <N-] (219) 


where now the "1" in (219) is again a continuous-envelope discrete 
mumcwon holding for all _n, o =n = N-1.§ In’general, in applying the 
method, since continued differentiations of a discrete function tend to 
introduce more and more discontinuities in the higher order derivatives, 
adjustments such as performed in (219) greatly reduce the "book-keeping" 
required. 

If now (219) is differentiated a second time the result is simply 

g‘‘(n)= o-N65 yy) oes ns N=! (220) 


Using (206) and (207) gives for (220) in terms of FFT expansions 


N-1 abe 2 ie. n 
> G(k) { 1-e } : 
k=o 
l N-1l a a fect 
=— y |-N(1-e ) Je (221) 
N 
k=o 


(221) must be true for arbitrary n which requires that the coefficients of 


exp (27 kn/N) vanish for all k#o0; hence 





. 2a aa me 
G(k) = ° ay eee oO eee N=-1 . (222) 
jie | 
The coefficient for k = o cannot be obtained from (221) since the 
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bracketed term vanishes; this is to be expected since k=o represents the 
constant term of the expansion and this must vanish upon differentiation, 
The k=o term can be however, trivially obtained as the mean value of the 
function over the interval, 
i.e., from (194) with k=o, 
-l 


N 
Dacia. (223) 
n=o 


Glo) = 4 


The results given in (222) and (223) have been checked using several 
Pre 29] programs to at least six decimal places, 
The procedure for the use of the method is now clear; if the 
discrete analytic-envelope function is polynomial in form it is evident 
that continued application of (202) must eventually result in the vanishing 
of the discrete analytic-envelope portion of the function leaving only a 
series of Kronecker deltas of various orders whose FFT coefficients can 
be written down by inspection. 
2. Coefficient of Recursive Waveforms -- The Exponential Function 

The question of how to treat infinite expansion polynomial 
analytic-envelope functions is also of interest. This would include func- 
tions such as the sine, exponential, etc. 

As an example of the the treatment of these functions consider the 
problem of finding G(k) when 


g(n) =eus ; o sn s N-] (224) 


g(n) being periodic (see Fig. 7) 


g(n +N) = g{n). (225) 
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Applying the derivative theorem gives 


aa) = ana __7 (n-lia 


_.7 (N-la 


aa ) 6, (n) ' a Oo (226) 


This can be expressed for all _n in the interval as 


g'(n) =o Pe NA 4 Qe NN) gn) - (1-e%) 8 (n) 


N-l)a 


=e M10) + (2 NN) gin), (227) 


or 
g‘(n) = g(n) (1-e°) “ sien) 6, (n) PaO ar th: = N=) (228) 


(228) is typical of the behavior of non-finite expansion polynomials for 
continuous functions: the derivative of such functions can be written in 
terms of the original function or certain linear combinations of the previous 


[28] 


derivatives The analogous result can be verified to hold true for 
discrete functions to first order; i.e., since the derivative in (202) 
includes only linear difference terms and ignores in (203) contributions 
to g(n At) proportional to ( At)” mm ( At)? , etc., the general expansion for 
the m-th derivative g ™ in) 


gm) gf) 


(x) =h ofa) + hom) +...b in), (229) 


] m-1 


where ho , eae hi; , are constants will deviate seriously from the 


1 
analogous expansion in the continuous case unles At is small. However, 
though the form of the expansion (229) is different from that for continuous 


functions the result (229) is completely correct for discrete functions if 


the constant terms are correctly taken into account. Such functions are 
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thus linearly recursive. If (228) is now expressed in terms of the FFT 


the result is 











2 2 
3 { (1-e ye (lee, ) } G(k) e 
k=o | 
27k 
N-1l | aa. a0 
i = 
=— { e*(1-e } aan ‘ ‘(230) 
N 
k=o 
This requires that 
a -Na) 
AME (l-e 
G{k) 2m7k 
(caer) 
i-e Na 
- as eS. ae a K=1 ,2,3,..33e N=] ° (231) 
ete 
Nieto N )y 


This particular result can be checked exactly analytically, for the function 


(224), when substituted into (194) becomes 





2 
N-1 -j n 
> ei na. N 
n=o 


a\r 


G(k) 





Zi 


n=O 


Ne 1 . 
(232) 


i 
N n=O 


where 





27k 
C6 (233) 


Equation (232) is simply a finite geometric series with the sum 
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(N-1)+1 
l ho5% 
G(k) = N (1 -x) 


N 
] -x 
ae ) 


l-x 


] 
DW CS 


1- -Na 
eo a NR) (234) 


)) 





ae 21k 


N(l-e N 


Note that in this case the term k=o has been obtained along with the usual 


other terms: this is because in (22) the derivative has been expressed 


-na 
recursively in terms of the complete original function e ; 


F, GENERAL PROCEDURE 

The above two examples illustrate the general procedure to be followed 
for finding the coefficients of a Fast Fourier Transform, The procedure 
may be summarized and generalized as follows: 


(A) If the analytic-envelope is finite polynomial in form contin- 
ually differentiate the function until only Kronecker delta 
functions remain. After each differentiation, adjust the 
analytic form so that the particular derivative holds at all 
points within the interval by extending the analytic-envelope 
over the entire domain o <n <s N-1 with suitable modifica- 
tion of the Kronecker delta amplitude. 


(B) Equate the spectral representations for g(n) to obtain the coeffi- 
cients for all k #oterms,. The k=o term is obtained as the 
average value of the discrete waveform. This can be obtained 
in many cases exactly by use of the finite arithmetic or geo- ~‘ 
metric sum formulas. 


(C) If the analytic-envelope consists of sub-regions over the inter- 
val o sn s N-1 defined by two or more different functions 
gj(n), go(n), etc., the problem can be handled by multiplying 
the appropriate sub-region function definitions by unit step 
functions for the appropriate regions. Differentiation of the 
resulting products using the usual product rule will lead to a 
simplified analytic-envelope function plus the necessary 
Kronecker deltas, The necessary work quickly becomes pro- 
hibitive unless the simplified analytic-envelope can be easily 
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extended over the whole interval as done in the previous 
problems. 


(D) If the function is transcendental, use the recursive procedure 
shown in the second example. 


In conclusion, the method presented in this section provides an 
analytic method for determining closed form formulas for FFT coefficients. 
These expressions can be used in theoretical discussions to study the 
behavior of the coefficients for various sampling intervals, size of 


discontinuity present, accuracy of computer programs, etc. 
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VI, TWO APPROACHES TO THE OPTIMUM ACOUSTIC ARRAY PROCESSOR 


A, INTRODUCTION 

A subject of great importance and deserving of continued study is 
the detection of acoustic communication signals in the presence of a 
background noise tielge? 9! In this chapter several approaches to the 
solution of this problem are reviewed and the results of simulating one of 
the approaches on a computer are reported and discussed, 

The usual means of detecting and receiving acoustic signals, 
corrupted by noise is by a spacially-distributed array of acoustic sensors, 
Such arrays serve the fundamental dual purposes of providing multiplication 
of individual sensor gain and spacial correlation. A typical array of 


sensors together with the associated signal processing network is shown 


moerig. 8. 





Array 






x, (t) 






Processing 





Network 
A 
s(t) = ? 


x, (t) 


x, (t 


Fig. 8. Array with associated processing network, 


Although analog, digital, or a combination of techniques may be invol- 


ved in the signal processing network in Fig. ,8, a few qualitative 


87 


generalizations can be made. If the input signal s(t) is wideband and 
quite high information rates are involved the decision to "go analog” is 
almost certain. The required computation rates for real time processing 
are very large in such cases even taking into account the possibility of 
using such tools as the FFT algorithm;thus digital processing is restricted 
to relatively low frequencies. However, the FFT algorithm, coupled with 
very high speed special purpose computers, have made possible data rates 
of 10 kHz and higher so that in the near future the use of digital real-time 
parallel processing may be expected to increase consideuatle ie 
There are many different approaches to the problem of optimum array 
processing; two main approaches will be discussed here. The first approach 
discusses the possibility of modifying classical analog optimum quadratic- 
cost-function detection and estimation theory to make use of the FFT 


algorithm; the second approach considers the use of Widrow adaptive array 


processor and possible modifications. 


B, THE MODIFIED OPTIMUM DETECTION - ESTIMATION ARRAY PROCESSOR 
[32 ] 


Following Cox closely let it be assumed that the L sensors in 


Fig. 8 have inputs x, (t), i=1,L and that the i-th input can be written as 


x(t) = [% m(t-r)s(r) drt+nt) , i=l,2,.0..,L (235) 
1 1 


-o© 
where m, (t) is a known linear transformation which takes into account the 
i-th time delay of arrival of the scalar signal s(t). m, (t) thus provides 
the generality to account for local dispersive properties of the medium, 


non-uniform array spacing, and phase distortion of particular sensors, The 
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function n (t) is assumed to be additive medium noise. The relation (235) 


can be re-written as the vector equation 


x(t) = | m(t-7)s() dr+ nit) (236) 


-© 
where now X, m, and n and Ir-dimensional column vectors. Several inves- 
tigators have formed large hyper-vector array ensembles of the quantities 
x, m, and n in the form of correlation matrices at various times the toe 
etc., in order to then use optimum estimation techniques 30,31], these 
procedures without the assumption of stationarity are confined to very low 
frequencies because of the great magntiude of the computation 
requirements. ee 


The assumption of stationarity will be made in the following sections 


in order to apply frequency domain techniques. Let it be assumed that 


mn) = (5) ues at mest Jo en” it 

=-—©o k=-o 
, at)=5  nik)el?™ *4yt (237) 

k=-o@ 
m(t) = - (fe at , (238) 

where 

ae 239 
Pi ase, (239) 


Then, substituting in (236) requires that 


O=E  {W'tk) - nf tk) eM F  yiga’ Kel? 


k=-o 


| = dt et T ee 


ann OO 
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=> {('k) -n'&) os er (Eo! (k)ed Ht 
k=l]2 —© 


6 (f-£, k) af } 


ae (2.40) 


ll 
7, &4 8 


{ ®‘(k) ~ wk) 0” (k) - nk) 


oo 
This can vanish identically for arbitrary t if and only if the bracket 
term vanishes which gives 

wv’ (k) = p' (kf, ) a’ (k) + n° (k) : Sidhe (241) 


(241) suggest arranging the equation in the matrix forms 


wy’ (1) n(1) a’ (1) 
g=(/ 22) gp a(2@), gale @ (2.42) 
BH (f,) © ° ns 
_[o HB’ (2f,) o 
E* ; ; w(3t,) . 
= diage( & (f,), a (En), 2.) (243) 
Or 
b=potyn (244) 


where wv’ (1) represents the Fourier series coefficients of the entire array 


of L elements at frequency fie w‘(2) the Fourier coefficients of the 


entire array at frequency 2f,, etc. If the signal s(t), input vector x(t) 


l a 
and noise n(t) are not band-limited the number of rows of the equation 


(244) is infinite; the usual assumption is that x(t) is band-limited to some 


er fr e 
upper frequency Lg where 
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f 
max max 1 


K max 
eee 249 
; (2.45) 


i 
~ 


and T is the observation period. 
Denoting the noise auto-spectrum at frequency kf, by R(k) where 
Rk) =E {ntk) nk) } , (264) 
and assuming that the period of observation T is large compared to the 
maximum lag for which the auto-correlation function is significantly 


different from zero, the cross spectrum R(k, satisfies 


E {n (k) n* (4) 


II 


Rik, £) 
=o a es (247) 


the noise covariance matrix is thus of the form 


RQl) 6 OQ «+ » 
On RIZ) ao 
-* ° eo R(3) 
= diag (R(l),; RZ), RG) a... (2 48) 


In (246) n* is the complex conjugate transpose of » . The matrix R 
is square and each submatrix on the diagonal is of dimension LxL. 
Suppose s(t) is assumed random with zero mean; a spectral covariance 
matrix P where 
P= Efotk)jox(g} , (249) 
can be defined. Under the assumption again that the expansion period T 


is long, P is also diagonal having the form 


oil 


] 
P = O Po O 
eee 
= diag (PD) sPysPaeee) a (250) 


where Py is the power in the spectral line kf Because of diagonal forms 


l e 
of the matrices yw, P, andR, the test statistic for the detection problem, 


C ,which when compared with a threshold statistic determines whether a 


signal has occured or not, can be expressed in the simple form 


] K / -] / 2 
C=5° 5 talk) (Hh R)* RE) Bw, , (251) 
k=] 
where 
p 
| ack) | = —4A——___-__.__—_, (252) 


1+p ( w’ (k) )*(RQK) B'(K) 


The form (251) is simply a sum of terms over the spectrum of the 
signal and is derived as the ratio of the conditional probabilities 


p(w / (worn)) 


C : 253 
p(w/ n) ae 
where 
p(v/(yotn)) = probability that during an observation interval 
T the input spectra contains the signal plus 
noise spectra, 
p(b/n) = probability that during an observation inter- 


val T the input spectra contains only 
noise spectra, 


If now a vector spectral function 9*(k), 
@* (k) = w'*(k) (R(k))™ (254) 


g2 


and a scalar spectral response function Q(k), 
Q(k) = a(k) @*(k) Ww (k), (255) 


are defined Parseval's theorem gives the simple result for C (T) 


k 
1 T | max 
C (T) = OT | | 2 Q(k)e 


Oo ~k 
max 


Jaet | oat (256) 
Equation (256) in practice can be accomplished by averaging the square 
law detected output of the spectral function Q(k). The function 6*(k) 
performs the operations of pre-whitening the input signal and then matching 
to the assumed known spacial structure of the signa. Fig. 9 shows 

the complete system; it is composed of a system of weighted bandpass 
filters on each sensor followed by spacial spectral summation, multipli- 
cation, square law detection and averaging. 

This processor can be generalized to supply with but slight modification 
not only detection and optimum estimation of a random Signal s(t) but also 
detection and optimal estimation for the known and unknown signal 
problems, Figure 10 shows the general all-purpose array processor assuming 
analog processing. 

The use of an analog spectrum processor has several disadvantages 
even though high real-time data rates are possible. Particularly objection- 
able is the fact that the banks of narrow band-pass "comb" filters required 
at the outputs of each sensor are difficult to maintain in precise adjustment 
track between sensors and are quite difficult to shift in frequency. In 


addition the process of analog multiplication is rather crude, 


I3 


As an alternative to the analog spectrum processor, the possibility of 
using the discrete FFT algorithm to perform the spectral analysis was 
considered. The essential modifications in the above analysis are as 
follows. 

Consider equations (237,238). Assuming stationarity the quantities 
x(t), n(t) and s(t) could be discretized by a digital sampler at the output 
of each sensor to provide the data values x(i At), n(i At), s(i At) where 
o <i s< (N-1) and the sampling interval At in any observation interval 
T satisfies 

N Ast =aeT 
Since the maximum harmonic that can be represented by the finite discrete 
N point Fourier transform is 


a N 
ae aaa fs 


io 
C= 


IO 





Te a 
the output of each sensor must be low-pass filtered to exclude all com- 
ponents of frequency greater than uae if the finite band-width of the 


sensor has not already provided the needed filtering. Failure to include 


such filtering will cause aliasing, 26] The discrete equivalent of (237) 
is thus 
N-1 jAEk ¢ N-1 | | oa 
x(€)==Z ke »n(f)=2 xn (ke ; 
k=o k=o 
N-1 — é 
( &)e =a wc (kK)e , § = 0,l,...., N-1- ay 
k=o 
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In (238) the spacial weighting vector m(t) is represented by a bi-lateral 
infinite Fourier transform and not a Fourier series. Ignoring this important 
difference for the moment and assuming that m{( &) can be represented 


by the finite discrete transform 


eee 


M N 
m(é)=ZD yp (gJe c= oO, , ee | (258) 
£ 


the discrete equation equivalent to (240) becomes 








N-1 jATE ewan = jen 
o=E {(u'k)-n'&))e 5 5 y'(de 
k=o T=0 Oo 
oe - 
ain } 
=E {(v'k)-n' We -E p'(ae 
k=o L =o 
N-1 j=27 (k-2) 
Calloies  e } 
T=0 
Zaks 27 2 
N-1 j € N-1 j 
=r {(w'W-n'wre % “-5 ye 
ao £ =o 


a’ (k) 6, (k-2) } 
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ay as g 


{ (k) - p(k) 0° (k) - 9 k) be * 4 (259) 


Il 
4 iS 


O 
Employing the same argument as in (240) the result 
w(k) = p(k) o (k) + n/(k) rakzo. ae peas (260) 

follows. Thus if (258) is valid the problem is again cast into the framework 
of (241) and the subsequent analysis follows directly. 

Now equation (238) in the analog formulation must involve a continuous 
frequency spectrum in order that any possible phase shift corresponding 
to any possible time delay 6 (t-t, ) relating the response time at the i-th 
sensor to that at some Convenient origin may be represented. The 
Kronecker delta in (259) does not possess such resolution; it can represent 
at most discrete time delays, separated by a minimum time of At sec, and 
a maximum time of N At=T sec. Thus only certain discrete angular resolu- 
tion is possible with the L element array under the representation restric- 
tion (238); the array has an angular boresight error of 


C At 
6o~T5X 


(261) 
where c is the local acoustic velocity and 6X is the fraction of a 
wavelength between elements for a uniform array detecting acoustic signals 
of wavelength X . Hence for fine angular resolution At should be small. 
A more fundamental restriction is imposed by the maximum time delay 
representable i.e,, NAt. Typical high resolution acoustic arrays may be 
of 30 Ato 50 X in aperture and in an endfire configuration this means that 


time delays of the order of 


delay ~ a (262) 
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across the array are possible, As a consequence even if T were not required 
to be long in order to guarantee the diagonal properties of the matrices P 
and R, it would be necessary for it to be long in order to satisfy (257). 

The net result of these considerations is that provided T is at least 
greater than some criterion such as (262) and At is small enough to satisfy 
desired angular resolution requirements such as (261) the analysis in (259) 
is valid and hence the FFT digital spectral analysis can be applied to 
optimally detect and estimate in the least squares sense the presence of a 


signal s(t) in the presence of background noise. 


C,. ADAPTIVE ARRAYS 
1. The Widrow Scheme 

The approach outlined above is classical in nature and relies on 
the statistical properties of signals alone to perform detection and esti- 
mation. The procedure requires the inversion of the noise covariance 
matrix which involves considerable computation and can only be com- 
pletely avoided if the noise is known a priori to be white. Further, apart 
from increased gain (which in general amplifies both the signal and noise) 
no particular correlation properties of the array are utilized. 

Conceptually, as previously stated in the introduction , an ideal 
array processor should be able to adapt to both the specific spacial- 
location medium properties where the array is to be placed and an incoming 
signal, 


[40 ] 


Widrow and his co-workers have devised an adaptive scheme 


which appears to have some of these properties as wellas other advantages 
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over the classical procedure in those cases in which the statistics of the 
signal s(t) are known and the noise field has non-uniform directional 
properties. No particular properties of the noise field are required except 
that the incident noise direction be different from that of the desired 
incident signal. The initial development given below is considerably 
expanded over that included in Widrow's paper, being a generalization 

to the vector case of a paper by Mantey. 41] 

In Fig.ll is showna discrete sampled-data network involving 
time delay networks for the i-th array element shown in Fig. 8. From 
Fig. ll, the output Y, (k) can be written as a function of the previous 
inputs and outputs as 

Ms (k) = ary, (k-1) + ay (k-2) tees tay, (k-N) 


+b x(k) + b.x,(k-1) +.... +b,,x.(k-M) , (263) 
oi li M i 


where x, (k) is the input sample at time kAt, x, (k-1) is the input sample 
at time (k-1) At, va (k) is the output of the network at time k At, y, (k-1) 
is the output of the network at time (k-1) At, etc. The coefficients a, 
oee-,a,. are referred to as feedback coefficients and the coefficients 


N 


Do. rar Di, are called feedforward coefficients. The term feedback here 
comes from the fact that the output i (k) is determined not only by the 
present and past inputs but also by the past outputs of the network, Such 
networks thus possess "memory" in that, indirectly, the output Y, (k) 
depends on all previous inputs x, (k) of the network. This is a highly 


desirable property. In detection theory the use of the matched filters to 


estimate threshold levels optimally in the least mean squares sense is 
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[32 ] 


important * such a filter is required to have an impulse response which 
is the convolution of the known signal waveform which is to be detected. 
Thus if the network in Fig. 11 is required to simulate such a filter itis 
necessary that its impulse response be as close to the waveform assumed 
a priori to be contained in x, (t) as possible. An important consideration 
in this regard is the relative efficiency in providing memory of each of the 
terms in (263); the impulse response is non-vanishing for the purely 


feedforward terms bos b, peeee,D.,, Only over the intervals (o <kAt < M&S) 


M 
whereas the feedback coefficients provide semi-infinite non-vanishing 
response, (o skAt < ~), Thus a much simpler feedback network can 
provide theoretically a sampled impulse response which can be made 
more nearly a "match” for the desired signal. As pointed out by 


[57] 


Mantey the cost in computing time and not computer storage is the 
essential limitation on the real time application of digital signal proces- 
sing schemes; it can be seen from (263) that each feedback section "costs" 
the same as each feed-forward section in computing time. The sample 
data network containing a single input b x(k), and N feedback terms 

any ¥(k-N), ao -a,¥(k) is called an auto-regressive network and has been 
extensively studied by Whittle 42 | and Claerbout 43] who conclude that 
the network has apparent advantages over other schemes in that it appears 
to yield more physically realizable models for stochastic processes which 
are time variable. 


Referring to Fig. 11 and equation(263) it is clear that the output of the 


entire L element array can be written in vector form as 
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N M 
y(k) =2 A_y(k-n) +% B_x(k-m) (264) 
n=! m=o 
where 
Y, (k) x, (k) 
yk) = : ; x(k) = ° (265) 
y, (k) x, (k) 
A, = diag (1 8on! ves aye B = diag(b, sb, os — 1b, : 
(266) 


In (264) there are N LxL matrices A, and M+l LxL matrices Be 
The theory of sample-data systems may be most elegantly treated 
[44] 


by the Z transform , which is defined in terms of a complex variable 


Z anda uniform sample sequence f(k) as 


F(z)=Z_ £f(k) ae 
k=o 


= 2(f(k)) (267) 


and has the property that the transform of a unit pulse delayed 2samples 


in time Oy yg , is 
ee 
Z( By? =z”. (268) 
Applying the transform to the system (264) thus gives 
N M 
Y(z) = 2 AT 4 (y(k-n) +2 Bj 2&(k-m)) 
n=] i= 
a -n _ -m 
= owl. 2  Wiz)i+ see Bez X(z) 
n = m = 
n=] m=o 
" -n -m 
= (Az) Y(z)+ 2 (BZ ) Xz), (269) 
n=! a m=o 
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where 
¥(z) = B(y(k))  , =X(z) = 2k) . (270) 
Defining the matrices A and B by 
N 


A=D Az 
n 
n=l 


n 


N 
= 2 diag(a,_.a 
n=! 


jz 


Nene oad 6 


N 
-n =a =n. 


N _, N 
= diag (S a Zz Beers a Z eae eS a 2 ), 
In 


M 
B=Z diag(b 
m=o 


ieee aan cae ee Ta 


M 7, OM _ M 
= diag (ZX bet BZ eee DB bz), 
m=o m=o ii=—© 


the equation (269) becomes 
Y¥(z) = AY(z) + B X({z), (272) 
which may be written as 


y(z) = (-a)7! B X(z) 


-m 
N 
l+Z a_z 


ie ln 


M 
ub? 
m 


= diag( op ete eONay (273) 


—=T}: 
The transfer function for the i-th sensor is thus 
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-1 -2 -M 
ae b +b, ,z +b) 52 teoe tb ,,2 


X, (z) = -2- (274) 


l-a aan Z se 
Ea 12 ce 





G (z) = 
i 

1N 

The denominator polynomial in (274) determines the stability of the network; 


[44] 


it is shown in Ragazzini and Franklin that if the magnitude of the 
complex roots of the denominator polynomial are less than unity the net- 
work characterized by equation (263) is stable. Since the transfer 
function corresponding to (273) is diagonal and involves no cross-coupling 
between the individual sensor networks, the question of array stability 
depends only on the stability of each individual sensor; if all of the 
sensor networks individually are stable and have finite bounded outputs, 
their sum will also have finite bound and be stable. However, the 
problem of stability in (274) is not necessary to consider at the stage of - 
the problem develpment here since the general sensor network in Fig. 10 
will be considerably modified in the subsequent development. 

The basic Widrow adaptive network is obtained from (274) by 


setting all of the feedback coefficients a i=l, .«i., N to zene, ive 


l a 


the network is a pure feed-forward network such as shown in Fig. 12. 
Fig. 12 represents a general M-th order feed-forward network and 
Widrow indicates that the identical delay networks shown will alow the 
network to receive signals over a band of frequencies if the weights 


Dye eoee, D.. are properly adjusted. In Fig. 13 is shown a single stage 


M 


of the network with the delay adjusted to 7/2 radians for some input 


frequency f Since at frequency f, the general output would be 


i it 
y(t) = b. cos(2mf,t + @) “ b sin(2mf,t + ¢,) (275) 
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it is clear that such a section could represent perfectly an input signal at 
frequency i, with arbitrary incident phase angle. Returning to Fig. 12, 
it is clear that M such sections in cascade will allow the exact repre- 
sentation of any input signal which can be written as the finite bandwidth 
Fourier series 


2 Oe (276) 


which represents a real bandwidth of Mf, Hz. If the delay is not pre- 
cisely m/2 radians for the lowest frequency of interest in the input 
signal the array can still compensate by readjustment of the weights 


a but the magnitudes of the higher order coefficients quite 


process Any 
rapidly (in the case of many delays) become very large. More will 
be said relative to this point later. 

Since the matrix B in (273) is diagonal it is convenient in 
considering the Widrow formulation to define a column weight vector w 


whose components are the diagonal components for the i-th array sensor 


from B before the application of the Z transform, 


aoe (277) 


Although the whole array of sensors could be handled corporately, it 
is only necessary to consider the sensors individually since there is no 
cross-coupling between sensors during the adaptive process to be 


described, 
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Defining a vector x, (t) as the array of inputs to the weighting 
coefficients in Fig. 12, the output Y, (t) can be written as 


y, (t) - w*(t) x, (t), (279) 


or in discrete form 


y. (k) = w*(k)x,(k) . (280) 
1 aL 


A very necessary requirement for the use of Widrow's adaptive algorithm 
is knowledge of the signal s(t), or discretely s(k). If s(k) is known an 
error signal e (k), can be defined as the difference between the value 
s(k) and the present output of the network given by (280), 


€ (k) = s(k) - w/(k)*x, (k). (281) 


A local gradient search technique, based upon the method of steepest 
descent, is then used to establish an algorithm for iterating the weight 
vector w(k) to a new value w(k + 1); the algorithm 


w(k+1) = w(k) - 2c) € (k) x, (k) (282) 


can be shown to converge after a large number of iterations to the Wiener- 
Holf solution. In (282) om is a proportionality constant which is always 
negative and usually satisfies 
Ie 
-l <c a x, *X, =O (283) 
i=) 

The necessity of knowing the signal s(k) which is to be received 
is the greatest restriction of the algorithm; since the signalis notin 
general known, except for possibly its second order statistics an a priori 
assumed artificial signal s(k) must be continually injected into the 
processor in order for it to adapt, This can be done in two ways. The 
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artificial signal s(k) can be injected into the processor directly for an 
adaptive period io during which it cannot be used for receiving, and 
then the artificial signal is turned off and the performance of the pro- 
cessor "decays" to degraded performance over a period oe during 
which the processor can be used for reception. The mechanism for 
doing this is shown in Fig. 14. If the power spectral properties of the 
artificially injected signal s(k) closely approximate those in the incoming 
waveform xX , &K) the of, periods will be long relative to the Tig periods. 
An alternative method is to provide two processors in parallel; an 
adaptive processor into which an artificial signal is continually fed to 
produce the optimum weights w to which a second processor is slaved 
to have the same weights and can thus continually receive external 
signals. Neither approach appears to be very satisfactory; the on-off 
approach is clearly useable only part of the time and the continuous, 
two-processor approach is expensive to implement and is "prejudiced" in 
favor of the apriori assumed signal s(k). 

2. A Visual Display Program 

Modifying a program, on which substantial previous programming 

effort had been contributed by LCDR William Moorehead of the U. S. 
Navy, an interactive study of the Widrow algorithm was performed, The 
program employs an IBM 360/67 computer with attached Model (2250-1) 
display unit; the flow chart for the program is shown in Fig. 15. By 
means of attention keys located at the display console (see Fig. 16) the 
program simulation could be interrupted dynamically at any point to 


provide additional information, change parameters, or produce permanent 
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Calcomp plots of any of the console displayed results. 

The basic information required for input to the program is the 
assumed direction and spectrum (center frequency and bandwidth) of 
the desired signal. The desired input signal power is automatically 
adjusted by the program to be unity with uniform power density being 
spread among the various sideband components over the desired 
bandwidth. The actual signal is then synthesized by use of the inverse 
discrete FFT algorithm. 

The noise field in the experiments that were performed was 
assumed to be incident normally to the array (see Fig. 17; as such the 
effective aperture of the noise signal is greater than that of incident 
signals for all incident signal directions not also normal to the array. 
Widrow employed sinusoidal "noise" for his test examples; the present 
investigation used random noise. The additive random Gaussian noise 
was added to the array from a random number generator program which 
produced random numbers of zero mean and known standard deviation. 

The time delay taps on the delay line were variable from a 
minimum of every sampling interval to any larger integral multiple of 
the basic sampling interval, 

At any given time in the adaption procedure up to four plots of 
directivity pattern vs. frequency and four plots of frequency response 
vs. angle could be obtained. This allows the receiving patterns at 
frequencies different from the desired frequency and frequency response 


diagrams for angles different from the desired look angle to be obtained 
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simultaneously. In addition the program provides the output spectrum of 
the array filter and the array filter weights w(k). 

In Appendix F are shown the results of running the program for 
the physical situation illustrated in Fig. 17. 

In general the results are as follows. The filter provides excellent 
Suppression for all frequencies outside of the band 400 <f < 600 cps. in 
the desired look direction of @= 45°. The response of the filter to other 
frequencies and angles of incidence is much reduced with the exception 
of frequencies above 1000 He where the pattern exhibits the familiar 
structure of linear arrays with element spacings greater than )A/2. i.e., 
the many lobe problem. This problem can be easily overcome, however, 
by either closer array element spacing or by ignoring the pattern and 
placing low pass filters of pass-band, o spass-band sat) where fy 
is the desired signal frequency of interest. 

In addition it is clear that the adaption time for the filter to 
reach steady state is quite long for general white noise inputs of the type 
used here; much longer than for the Narrow band simusSoidal noise 


[40] 


used in Widrow's paper. 
3. Other Approaches: Future Research 
The results reported in this paper relative to the Widrow filter 
Simulation program are minimal; rather it is thought more appropriate to 
make a few generalizations and suggestions for future research based upon 
results already observed. 
The Widrow filter must have knowledge of the signal to be detected 


(at least the nominal frequency and bandwidth); since a priori the signal 
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to be detected is of course not known an assumed signal s(k) must be 
applied to adapt the filter. From a practical point of view this must 
mean a continual process of adaption using an ensemble of maximum 
likelihood input signals s(k) as the adaptive signal inputs; the speed 
of adaption is thus of fundamental importance as all angles, band- 
width frequencies, and nominal center frequency ranges of the ensemble 
must be scanned during the acquisition phase of operation. The single- 
mode two-processor adaption algorithm is thus almost a necessity. 

secondly, the real time speed of adaption is strongly a 
function of the amount of adaption computations the processor is required 
to do. In this regard the alternative procedure of assuming every input 
is narrow-band to first order during the acquisition phase (and thus 


requires only the computations for Doe b, for any element) appears to 


l 
have special merit, especially for the detection of amplitude modulated 
signals, During the acquisition phase b, sod execene Diy would be maintained 
at zero, After acquisition of a signal the "fine structure" of the broad- 
band spectrum can be provided by adding other spectral components to 
the adaption signal and adapting the weights bib) bo, oe eee Diy as 
before. The signal output during the preliminary adaption phase would be 
quite distorted; this might be a tolerable price to pay for much faster 
speed of adaption. 

Another possibility for increased speed of adaption as well as 
providing more stable weights would appear to be the use of non-uniform 


time delays. There is nothing in the recursive algorithm (276) which 


requires that the delays all be equal; the algorithm requires only that the 


weight vector, error scalar, and particular input signal be known at any 
sequential time step. Consider for example the amplitude modulation 

case in which the desired signal is as shown in Fig. 18. This signal 

s(t) is wide-band in the sense that the (band-width) / (center-frequency) 
~.5 but is narrow-band in the sense of equation (276). If all of the delays 
are of size 1/4f, only the fundamental is very effective in simulating the 
input signal; the amounts of the higher order harmonics in (276) provide 
very small envelope corrections to the waveform synthesized by dD. and 


b This suggests the network shown in Fig. 19; the adaption procedure 


i e 
consists of first doing a first-order narrow-band adaption as previously 


described to obtain bo and b._ for the first delay of i/4i, with by mie a, 


i 


Dis fixed at zero and the desired adaption signal being a monochromatic 


signal of frequency f When the gradient of the error function €e (k) 


l e 
decreases to a pre-assigned threshold value, the additional signal 


components are added to the desired adaption signal s(t) and the weights 


Ee, D-, b 


z 1 2 are all allowed to adapt to steady state value. If 


Pes «9D 


M 


the desired signal s(t) is of the form shown in Fig. 18 only a very few 


additional weights dD, 1 cece, DO, Should be required and the adaption time 


M 
should be very rapid. 

4. The Use of Feedback Networks: The Kalman Filter 
[41] 


The work of Mantey indicates that it may be possible to 
produce a better adaptive synthesis procedure using feedback networks. 
However, he shows that the sum-squared quadratic error function 


5 Le (k)? =5 (s Gem, (284) 
1 1 
k=o k=o 
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possesses a unique minimum with respect to adaption only if the feedback 
coefficients ay pn enetevons any are all fixed, He then shows that modified 
adaption procedure in which the desired signal is used directly as the 
adaption signal applied to the feedback coefficients (Fig. 20), leads to 
a quadratic performance criterion which has a unique minimum, This 
new adaption procedure, in the case that x(k) is white, has also been 
shown to be a stable ames 4 122) 
Thus it would appear that perhaps substantial improvements in 

adaption performance might be realized by use of time variable feedback 
networks employing one or more feedback loops. The Kalman eiltere46+47 3 
is such a filter; itis a least-square-estimation filter which provides 
optimal estimates of the state vector s(k) recursively. The problem may 
be formulated here as that of obtaining the best estimate s (k/k), of the 
state vector s(k) at time state k given k previous observations x(k) and 
the previous optimal estimate s(k/k-1) of the state vector. Itis assumed 
that s(k) satisfies the dynamic system relations 

s(k) = ®(k,k-1) s(k-1) + A(k,k-1)n(k) (285) 

x(k) = H({k)s(k) + v(k) (286) 
where: 

s(k) is the state vector (unknown) at time state k 

€(k,k-1) is the state transition matrix (known) 

A(k,k-1) is a known transition matrix 

H(k) is a known transformation matrix 

x(k) is the observation vector 


n(k), v(k) are sequences of independent Gaussian random vectors 
with Zero means and the known covariance matrices 
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E{ nik) n*(]=Q6, , Elv(k) vte)I=R8, , 


L 

ELatkyy*(2)) = om (287) 
They represent assumed additive excitationand measurement noise signals 
respectively. 


[47] 


Following Lee, the optimal estimate s(k/k) is assumed to be 
the sum of the extrapolated estimate, based on past estimates anda 
weighting "gain" G(k), times the residual difference between the new 
measurement and the extrapolated previous measurement, 

8 (k/k) = ®8 (k/k-1) + G(k) [ x(k) - HS(k/k-1)] . (288) 
This can be re-written in the form 

ae oe ee) + Ge) (289) 
This suggests in the case of two equal time delays the state transition 
flow diagram shown in Fig. 21. The linear equations equivalent to 
(289) are, from Fig. 21, 


S, (k/k) = 14-9, 19,1 18, k-1/k-1) +14, 57 13, (k-1/k-1)4 9, x(k) 


311% 
(290) 


8 5 (k/k) = (41-9, 14 18, (k-1/k-L) +10, 9-9, 1) 9 185 k-1/k-1) 4g, x(k) 


(291) 
A -la 
$, (k-1/k-1) = z 8, (k/k) (292) 
A —-lA 
s,(k-1/k-1) =z "s, (k/k) (293) 
where 
S11 912 yo 11 Y%12 
G = | He , &= | 
So, 920 2 Go, %o/. 
(294) 
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Since the desired output is S(k/k), the system (290,293) can be solved 
for the transfer function 51 (k/k) / x(k), 


S, (k/k) - G1, U-(y-9,,4,)2 1+9,,[h,]2 —_ 


x(k) 7 ¥, 


= 
a ae (295) 


¥2 


where 


7 = =i 
H = 1-1-9149) V0-( 0-991 4 0)2 )- (9-9) 1%) 0) 


-2 
° (41-99 1% 12 


os -] 
Wo= UU 1-91 1H 2-Fa BoB — + EG) 1-94 111) 9-9 21 919) 
- (6,59, 6 )(d.-9,,H 12 ~ 
2 oe) Lae % ) 2" oh ; 
This, by comparison with (274) is equivalent to the sample data system 
shown in Fig. 22. For any desired signal s(k) it is possible to thus 
obtain the optimal least-squares estimate of S| (k/k) at the output of 


the network provided 


bo (k) = 9), (k) (298) 
b, (k) = Ld 9 (K)-, 1 (kK), 9 (k)-& 4 (kK) +9, , (kK) 4, kk) ] (299) 
a, (k) = Ld, 1 K-91 KA | kK) +0, (k)-9, | (kK) A fk) J (300) 


a, (k) = U9, | (K)-g, , (K) 4, , (K)) (@,, (k)-9., , kb, 5 kD) 
— (8 9 kK) 1 KO) 9 (KN) (8, | (K)-9, (kK), KD 
(301) 


If the signal process s(k) is stationary, the matrix © will be 


leh 


time invariant and the gain matrix G(k) will quite soon reach steady states; 
this leads to the exciting possibility that the steady state pre-calculated 


values for a,,a 


j bos b 


9! ; may be used directly in (295). Although this 
would not give optimal estimates during the usual transient phase of the 
adaption procedure, previous experience with the technique of using 
fixed gains indicates that it begins providing optimal estimates after a 


very short time. 


Dis 


Vi, CONCLUSIONS 

This research has concerned investigations into the inter-related 
topics of large-domain integration techniques for solution to the scalar 
multi- dimensional variable-parameter acoustic wave equation and 
adaptive antenna-array processors and algorithms. 

After formulation of the acoustic propagation medium equations to 
be integrated in properly posed form, their solution by means of various 
finite-difference integration techndiues was investigated. Among the 
results that are considered new or significant are the following. 


[45] 


A method for “tearing" or partitioning very large variable parameter 
domains in cases where steady state propagation characteristic are 

desired is given. Such positioning is essential if the integration of the 
wave equation for large domains is to be carried out on reasonable-sized 
machines or in background mode in a multi-tasking environment. Computer 
calculations verify that consistent results are obtainable using the 
method, but its practical use depends upon the efficiency of the integra- 
tion scheme chosen and "asking the right questions." Since, inherently, 
the method is capable of giving the complete steady-state pressure velocity 
distribution profiles at all points of space it is not surprising that quite 
large computing time requirements result. However, if the integration is 
confined by program control to only those regions of the spacial domain 
wherein the state variables are above an assigned threshold, great 
reductions in computing time are possible. Indeed there appears noa 


priori reason why ray-path analysis cannot be used for regions wherein 


the properties are known to be quite regular with a switch to direct 
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wave-equation integration when approaching regions of unpredictable 
behavior such as the incidence of convergence zones, highly variable 
parameter domains, defraction regions, or where three-dimensional 
behavior cannot be ignored, 

A new method for integrating the wave-equation in three dimensions is 
developed which uses the Fast Fourier Transform digital algorithm. The 
method is fast, accurate, and easily implemented for general rectangular 
domains but as with all Fourier techniques is restricted to constant 
parameter regions. The method is semi-implicit and appears to possess 
considerable stability even for fairly large integration step sizes, 

The Douglas and Gunn tri-diagonal integration scheme was imple- 
mented for the wave equation in three dimensions. The scheme is faster 
computationally than any of the other algorithms tried in the case of 
variable coefficient problems. However, it requires a great deal of 
programming effort, large amounts of work storage, and the conditions of 
stability are not precisely known. 

As an aid in understanding and applying the general discrete Fourier 
Transform, a new method for finding analytic closed form expressions for 
the coefficients is given. The method is basically an extension to discrete 
domains of the impulse function method of continuous distribution theory. 
Although the method does not guarantee that the investigators efforts will 
be small, it appears that coefficient expressions can be obtained quite 
easily for functions having relatively few discontinuities in analytic- 


envelope, 
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Relative to processors which seem especially appropriate for use in 
variable parameter environments the topic of adaptive array processors 
and algorithms was considered. First, the classical optimum statistical 
estimation procedure as it applies to array processors was reviewed, 

In the classical analog approach it is shown that discrete spectrum 
algorithms of the FFT type can be substituted provided in the analogous 
continuous spectrum variables proper precautions are observed relative 
to estimation and/or detection interval T and sampling interval At. This 
generally will require that N, the number of sampling points for dection 
interval T, be quite large; this is no particular restriction and in fact the 
FFT algorithm approaches maximum efficiency for large N. 

The second point of view is that leading to automatic adaptive 
synthesis procedures in which quadratic-error functional minimizations 
are carried out. The Widrow scheme is considered as a special case of 
a general feed-forward, feedback network of the type considered exten- 


sively by Mantey. Results of a visual display simulation program are 


reported and the use of other possible adaptive algorithms such as variable 


delays and the Kalman filter are discussed, 
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APPENDIX A 
THE DOUGLAS AND GUNN EQUATIONS 

The Douglas and Gunn equations iteratively solve the wave equation 
in three dimensions by computing a series of three approximations to 
the state components at the time (r+1) At in terms of state components 
at the time rAt. The particular equations to be solved for any one 
approximation are tri-diagonal in terms of the implicit variables; the 
solution to the implicit portion of the previous approximation acts as 
an explicit input to the next approximation. Thus the implicit system 
of equations to be solved at any stage are tri-diagonal and hence, in the 
case of the three dimensional wave equation the work of solving the 
system is only slightly greater than three times the work required to 
solve the equivalent one-dimensional wave equation. 

Denoting by u(Z ,m,n,r) the known state of the system at position 
L&x, mAy, n Az, r At, by u*(£,m,n,r+1) the first approximation to the 
solution of the wave equation at state (r+1) At, by u**(Z,mn,r+1) the 
second approximation, etc,, the equations to be solved are shown below, 
Por purposes of brevity only the indices that vary in any given term are 


explicitly written down, 


p*(r+1)-p= -@ lv* (4+1,r+1)-v* (f-1,r+1) + v (41) - v, (4-1) 
+2 (v (m+l) - v (m-1) +v_(m+l) - v (n-1)) J] (Al) 
y Y z Z 
vx (r+1)-v_ = -@ [p*(44+1 ,r+1)-p*(£-1,r+1l) + p(4+1) - p(&1) ] (A2) 
vr (+l) eas -2a@ [p(m+1) - p(m-1)] (A3) 
v *(r+1) a 2a {p(nt+1) - p(n-1) ] (A4) 


End of first approximation equations 


p**(r+1) -p = -a [lv *(£41,r+1)-v.*(£-1,r+1) + v_ (441) -v_ (4-1) 
+v**(m+l,rt+l) - v **(m-1,rt+1l) +v (mtl)-v (m-l 
y ( ) y ( ) a ) y! ) 


re 2(v_ (n+l) - v_(a-1) dadl (AS) 


v** (rt) Veet [p*(2+1,r+l) - p*(£-1,r+1) + p(£+1) -p(2-1) ] (A6) 
eer”) ro (p**(m+1 ,r+1)-p**(m-1,r+1) +p(£+1)-p(4-1) J] (A7) 
v ¥*(r+1) ke -2a [p(n+1)- p(n-1) ] (A8) 


End of second approximation equations 


= = — * -—ty xX a pe = 
p(r+l)-p aly? (£+1 ,r+1) id (2-1 ,r+1) + v, (£41) v,(é 1) 
+ ve prtl) - yy aed Tt avegim +L) -v Argel 
3 v_(m+1 iTt1)-v (n-1 T+1)+v, (nt1)-v, (n-1) ] (A9) 
Vv (+1)-v = -a (p* (£+1 ,r+1)-p* (4-1 ,r+1)+p(£4+1)-p(é-1) J (Al 0) 


Moe -a [p**(m+1 ,r+1)-p**(m-1,r+1)+p(m+l)-p(m-1) ] (Al1) 





Ve (r+1)-v = -a [p(nt+1 ,r+1)-p(n-1,r+1)+p(n+1)-p(n-1) ] (A12) 
where 
_ cAt : (Al3) 
Oo” UA AX 
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APPENDIX B 
THE SIGNAL INPUT TO A MATCHED LINE 
This appendix shows the response of a matched output one-dimensional 
domain to a step input function applied at time t=o sec. The domain is 


of unit length and is subdivided into 100 segments. 
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APPENDIX C 
THE COSINE INPUT TO A MATCHED LINE 
This appendix shows the response to the same domain as described 


in Appendix B but with a 1 Hz cosine signal input. 
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APPENDIX D 
TWO-DIMENSIONAL VARIABLE-VELOCITY PROBLEM SOLUTION 
This appendix shows the two-dimensional pressure distribution for a 
line source initial condition applied to a medium having a variable-velocity 
profile. The example shown is for a case in which a linear increasing 
velocity profile is present. This condition would be present in seawater 


as the depth increases, 
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APPENDIX E 
TWO-DIMENSIONAL NON-ZERO INITIAL CONDITION SOLUTION 
This appendix shows the transient solution to a deformed two-dimen- 


sional membrane problem. 
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APPENDIX F 


The graphical results of running a Widrow adaptive filter simulation 
program are shown in this appendix. The desired signal to which the 
array responded had the following characteristics: 


center frequency = $00 Hz 


band-width = 200 Hz and 0 Hz 
desired look angle = 45 degrees 
power Signal input = 1.0 watt 


The noise signal is of standard deviation 0.707 and is provided by the 
random number program GUASS; thus for wideband signals the total noise 
power input to the array may exceed that of the signal before filtering. 
The noise is incident at an angle of 0 degrees relative to the normal of 
the plane of the linear array. 

The filtering action occurs in both in frequency response selectivity 
and spacial angular selectivity. The graphs on the following pages are 
grouped into frequency response graphs and receiving pattern graphs. 
There are four frequency response graphs and four receiving patterns 
shown for each time of adaption. 

The frequency response graphs show the response as would be observed 
for the different angles of i a0 45°, 90°. The results show that the 
desired signal of 500 + 100 Hz is strongly attenuated for all directions 
other than the desired direction of 45°, 

The pattern plots show the spacial response of the filter at the 


frequencies 100, 200, 500, and 1000 Hz. The results show that these is 


greatly reduced array pattern gain for the desired frequency of 500 + 100 Hz 
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in all directions except the design direction of 45°, In addition, the 
response of array to the frequencies of 100 Hz and 200 Hz are greatly 
reduced in the design direction of 45°, 

The one prominent exception to the desired behavior is that the 
filter has appreciable gain and frequency response for signals or 
noise at 1000 Hz or higher. This is correctable as mentioned in Chap. 
VI by appropriate low-pass prefiltering. 

The results for zero bandwidth desired signals are very close to 
those shown for the wide-band signal except that the rate of adaption 


is considerably greater in the desired look direction. 
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Figy 1. Diffusion equation solution by Fourier series and 


stable discrete difference methods, €urves = 
Fourier solution; dots = déscrete solution, 
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Fig. 2. Instabilities in discrete diffusion equation 
solution, 
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I(x,t) = Re(I(x)e ) 
Fig. 3. The One-Dimensional Transmission Line Equations 
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Fig. 4. Geometric model of the integration domains 
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Fig. 5. Function having discontinuity at 
Nn =2 
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Fig. 6. Derivative of function g(n) =n 
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Fig. 7. Recursive example: The exponential 
function 
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Fig. 13. Narrow-band feed-forward delay network 
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Fig. 14. The on-off adaptive network scheme 
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Terminate Problem 
Fig. 15. Flow Diagram of Array Simulation Program . 
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Problem Data: 


Desired Signal Center Freq, = 500 Hz 

Bandwidth = 200 Hz 

Number of Elements = 1] 

Number of Delays/Element = 2 

Pattern Plots at Angles of 0°, 30°, 45°, 90° 

Frequency Response Plots at 100 Hz, 250 Hz, 500 Hz, 1000 Hz 


Fig. 17. Physical Layout of Array Simulation Program. 
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Fig. 18. Wide-band Amplitude Modulated Input Signal, 
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